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Abstract. In order to increase the efficiency of the computer simula- 
tion of biological molecules, it is very common to impose holonomic 
constraints on the fastest degrees of freedom; normally bond lengths, 
but also possibly bond angles. Since the maximum time step required 
for the stability of the dynamics is proportional to the shortest period 
associated with the motions of the system, constraining the fastest vi- 
brations allows to increase it and, assuming that the added numerical 
cost is not too high, also increase the overall efficiency of the simula- 
tion. However, as any other element that affects the physical model, 
the imposition of constraints must be assessed from the point of view 
of accuracy: both the dynamics and the equilibrium statistical mechan- 
ics are model-dependent, and they will be changed if constraints are 
used. In this review, we investigate the accuracy of constrained mod- 
els at the level of the equilibrium statistical mechanics distributions 
produced by the different dynamics. We carefully derive the canonical 
equilibrium distributions of both the constrained and unconstrained dy- 
namics, comparing the two of them by means of a "stiff" approximation 
to the latter. We do so both in the case of flexible and hard constraints, 
i.e., when the value of the constrained coordinates depends on the con- 
formation and when it is a constant number. We obtain the different 
correcting terms associated with the kinetic energy mass-metric tensor 
determinants, but also with the details of the potential energy in the 
vicinity of the constrained subspace (encoded in its first and second 
derivatives). This allows us to directly compare, at the conformational 
level, how the imposition of constraints changes the thermal equilib- 
rium of molecular systems with respect to the unconstrained case. We 
also provide an extensive review of the relevant literature, and we show 
that all models previously reported can be considered special cases of 
the most general treatments presented in this work. Finally, we nu- 
merically analyze a simple methanol molecule in order to illustrate the 
theoretical concepts in a practical case. 
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1 Introduction 

Since its first practical applications in the 1950's (see [T] and references therein for 
a detailed historical account), molecular dynamics (MD) has become a powerful and 
well-established tool for the study of a wide range of systems, including but not limited 
to condensed-matter materials j2] , fluids |3|4j , polymers and biological molecules [5] . 

Among the latter, one of the most important families of systems is that of pro- 
teins. Together with nucleic acids, proteins are one of the main responsible of the 
enormous complexity and versatility that living beings exhibit. Although they are 
mostly constructed as linear chains of only twenty different monomers (the twenty 
proteinogenic amino acids), and these monomers, in turn, are made of only five types 
of atoms (H, C, N, O and S), the finely tuned combination of these elements is capa- 
ble of producing remarkable 'nano-machines' that perform almost every task that is 
complex in biological organisms [BJ. Thus, it is not surprising that a great effort has 
been invested in the last years in devising new theoretical and computational methods 
to describe and simulate these important biomolecules. In particular, many groups 
are pushing forward very ambitious scientific and technological agendas to be able 
to perform longer and more accurate MD simulations of increasingly larger proteins 
by different approaches |7I8I9| . In this work, we have in mind the long-term goal of 
producing more efficient methods for simulating proteins, however, the formalism we 
introduce and the calculations we perform are directly applicable to any molecular 
system, as long as it size makes it manageable in present day computers. 

The efficiency of a MD simulation is customarily defined as a suitable relation 
(e.g., the quotient) between its accuracy and its computational cost [TO]. The overall 
accuracy can be more finely divided into that related to the physical model used 
to describe the system and the accuracy of the method used for implementing the 
dynamics. The computational cost, on the other hand, can be split into that associated 
with increasing the size of the system and the cost related to the time-propagation of 
the dynamics [TT]. The fact that all these issues are intricately coupled is easily seen 
if we consider examples such as the comparison between ab initio MD based on the 
ground-state Born-Oppenheimer approximation |lll2ll3ll4U5j and MD simulations 
using as energy functions the ones known as classical force fields |16ll7ll8ll9l20l21j . 
On the one hand, it is clear that the accuracy of the two physical models is not the 
same; on the other hand, the size of the systems that can be practically tackled is also 
different, since quantum chemical methods present a cost that scales typically faster 
with the number of atoms than that of force fields [22I23| . Another example of this 
interplay between accuracy and cost is the use of coarse-grained descriptions, which, 
by selecting as interaction centers larger entities than individual atoms, and thus 
changing the physical model, allow to reach larger systems and longer time-scales than 
atomistic simulations using either force fields or ab initio methods [24 25 26 27l28j . 

Yet another technique which can also be used to try to reduce the computational 
cost of MD simulations by modifying the physical model, and which is the object of 
this work (and of this special issue of The European Physical Journal), is the impo- 
sition of constraints on some of the degrees of freedom of the system. The concept 
of constraints itself is very general, and it can be used in many ways in the con- 
text of molecular modeling and simulation: For example, in Car-Parrinello molecular 
dynamics (29j . which is a type of ab initio MD devised to mimic ground-state Born- 
Oppenheimer MD, the time-dependent Khon-Sham orbitals need to be orthonormal; 
a requirement that can be fulfilled imposing constraints on their scalar product [30 . 
Also, we can use constraints to fix some slow, representative degrees of freedom of 
our molecular system, also called reaction coordinates, in order to produce free energy 
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profiles along them that would have taken an unfeasibly long time to be calculated 
should we have used an unconstrained simulation (see the Blue Moon Ensemble tech- 
nique in ref. [311 and also a number of related works in this special issue [32133] ). 
In this review, however, we are only concerned with the use of constraints in order 
to fix the fastest, hardest degrees of freedom of molecular systems (like, e.g., bond 
lengths or bond angles), with the objective of both reducing the effective size of the 
conformational space, and allowing for larger time-steps in MD simulations. 

The first objective is worth pursuing if we acknowledge that one of the charac- 
teristics of large, flexible molecules like proteins (and many others) is that they have 
an astronomically large conformational space; an issue that may hinder both our 
understanding of the problem and the possibility of designing algorithms that could 
efficiently sample this space 34 35 36J . These systems additionally exhibit movements 
with typical times that span a very wide time-scale, ranging from very fast bond vibra- 
tions in the range of tenths of femtoseconds, to the more interesting conformational 
changes related to biological function (like allosteric transitions or protein folding) 
which can take times of the order of the millisecond, or even the second [37138139] , 
Since the time-step in MD simulations has to be chosen at least one order of magni- 
tude smaller than the smallest typical time associated to the motion of the system if 
we want the results to be accurate |39I40I41I42I43| , this means that we need to calcu- 
late 10 13 -10 15 MD steps if we want to capture the mentioned biological phenomena 
in our simulation. 

The imposition of constraints that fix the fastest degrees of freedom holds the 
promise of alleviating both these two problems, however, as any modification of the 
physical model that describes our system, it may also carry with it a loss of accuracy 
that renders the simulation useless. In this work, we introduce the rigid and stiff 
constrained models in order to be able to quantify this loss of accuracy at the level 
of equilibrium statistical mechanics. In doing that, we benefit from the occasion to 
thoroughly review the relevant literature, and to try to unify the different, sometimes 
conflicting, vocabulary and mathematical definitions used to get a handle on the 
physical concepts. This endeavour is facilitated by the fact that the formalism used 
in this review is, as far as we are aware, the most general one in the literature, therefore 
allowing to obtain all previously discussed models as particular cases in which the 
(often implicit) approximations can be clearly identified. 

In sec. [2j we introduce the theoretical framework, beginning by the notation in 
sec. |2.1| and the Hamiltonian dynamics and statistical mechanics of unconstrained 



systems in sec. 2.2 Then, we describe and physically justify in sec. 2.3 what types 



of constraints can be applied, both in terms of the set of coordinates chosen to be 



constrained (sec. 2.3.1), and in terms of whether or not their equilibrium values are 



considered to depend on the conformation of the molecule (sec. 2.3.2); we term the 
constraints flexible if this dependence exists, and hard if the constrained coordinates 
are assumed to take constant values. In sec. |2.4[ we introduce the four possible com- 
binations of constrained statistical mechanics models used in this work and in the 
literature: the stiff and the rigid model, with either flexible or hard constraints. 



In sec. 2.5 we comment on the mechanism for comparing the different models and 
introduce the so-called Fixman potential to do so; also, we collect the different ap- 
proximations made in the literature with respect to the most general cases discussed 
in this work, and we establish the relevant relationships among them. In sec. [3J we 
provide a numerical example in the simple methanol molecule in order to illustrate 
the concepts, and we benefit from the occasion to also review the previous numerical 
analyses in the literature. Finally, in sec. |4j we outline the main conclusions of the 
work, and we mention possible lines of future research. 
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2 Theoretical framework and relation to previous works 

In this section, our objective, which is in fact one of the main objectives of this 
review, is to lay down the needed mathematical formalism to analyze the equilibrium 
statistical mechanics of constrained molecular models, and to do so in the most general 
terms possible. 

Of course, we are not the first ones discussing constraints in molecular modeling. 
As with any other scientific topic, there is a significant degree of arbitrariness in any 
attempt to fix its historical beginning, but we can talk about a discussion that has 
been going on at least for three or four decades at the moment of the writing of this 
manuscript. 

Therefore, what is the interest of a new account? 

We are confident to be able to convince the reader in the next subsections that, 
despite the many great works that have dealt with the problem so far, no article 
defines the problem with the same generality that we display here. This general 
formalism, which contains all previous descriptions as particular cases, allows us to 
clearly identify the different approximations that have been assumed in the literature, 
to relate the many models with one another, and to propose a consistent wording and 
notation which could facilitate the clarity of future developments. 

As we move along, we will indicate at each step which other choices have been 
made in previous works and how they relate to the general setting. 



2.1 Notation 

The system of interest (which we could call the molecule) is a set of n mass points 
termed atoms. The Euclidean coordinates of atom a in a frame of reference fixed in the 
laboratory are denoted by r a (Vectors of 3 components in the Euclidean 3-dimensional 
space are denoted by bold symbols.), and its mass by m a , with a = 1, . . . ,n. However, 
when no explicit mention to the atom index needs to be made, we will use r T := 
( rAt )^=i to denote the (row) TV-tuple of all the N := 3n Euclidean coordinates of the 
system. By r, we shall mean the corresponding column TV-tuple, or just the ordered 
set (r 11 )^^ when matrix operations are not involved. These coordinates parameterize 
the whole space, denoted by W, and the masses iV-tuple, m T :— (m^)^ =1 , in such a 
case, is formed by consecutive groups of three identical masses, corresponding to each 
of the atomaS 

Apart from the Euclidean coordinates, we shall also use a given set of curvilinear 
coordinates (also called sometimes general or generalized) , denoted by q T := (q^)^ = i, 
to describe the system. The transformation between the two sets and its inverse are 
respectively denoted by 



1 Some of the quantities appearing in the formalism, such as the velocities r M , the displace- 
ments <5r M , or the mass-metric tensor G^v (see the following sections), change like tensors of 
some type under a general change of coordinates; in such cases, the position of the indices 
is dictated by the tensor type (1-time contravariant tensors in the case of or 5r M , 2-times 
covariant in the case of G^ u ). Some other quantities, such as f M , F M , or r M , do not change 
like any type of tensorial object under a general change of coordinates; in such cases, the 
position of the indices is chosen according to notational convenience, or to the way in which 
they transform under some particular family of changes (e.g., linear ones). 



r" = R»(q), n = l,...,N, 
q» = Q»(r), fX = l,...,N, 



(la) 
(lb) 
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and we will assume that, for the points of interest, this is a proper change of coordi- 
nates, i.e., that the Jacobian matrix 

has non-zero determinant, being its inverse 

Although the curvilinear coordinates q are a priori general, it is very common to 
take profit from the fact that the typical potential energy functions of molecular sys- 
tems in absence of external fields do not depend on the overall position or orientation, 
i.e., they are invariant under overall translations and rotations, in order to choose a set 
of curvilinear coordinates split into q T = (e T ,u> T ), where the first six, e T := {e A )\ =1 
are any external coordinates that parameterize the aforementioned global position 
and orientation of the systenrl The remaining TV — 6 coordinates w T :— (w a )^ =7 are 
called internal coordinates and determine the positions of the atoms in some frame 
of reference fixed in the system. They parameterize what we shall call the internal 
subspace, denoted by I, and the coordinates e parameterize the external subspace, 
denoted by £ ; hence, we can split the whole space as W — £ x 1 (denoting by x the 
Cartesian product of sets). 

The reader is probably familiar with the use of frames of reference fixed in the 
system based on the center of mass and the principal axes of inertia, but, although 
this is well adapted to the parameterization of the dynamics of rigid bodies [H] , it is 
less so to flexible entities like molecular systems, where the relative distances of their 
constituents are not constant in time. Following ref. [45], we define a more suitable 
frame of reference fixed in the system to perform some of the calculations, although it 
is worth stressing that most of the mathematical formalism introduced in this work 
should not depend on this choice as long as the external and internal subspaces are 
well separated. 

To define this frame of reference, we select three atoms (denoted by 1, 2 and 3 
in Fig. [IJ in such a way that o, the position in the laboratory frame of reference of 
the origin of the frame of reference fixed in the system, is the Euclidean position of 
atom 1 (i.e., o := ri). The orientation of the frame of reference (x' , y' , z') fixed in the 
system is chosen such that atom 2 lies in the positive half of the z'-axis, and atom 
3 is contained in the (x', z')-plane, with projection on the positive half of the x'-axis 
(see Fig. [I]). The position of any given atom a in the new frame of reference fixed 
in the system is denoted by r' a , and we have that the previously unspecified external 
coordinates are now e T := (e A )\ =1 = (o x ,o y ,o z ,4>,9,ip), where (<f>,6,7p) are three 
'Euler angles' that parameterize the orientation of the 'primed' axes with respect to 
the 'unprimed' ones. 

More explicitly, if E((f>, 9, ip) is the Euler rotation matrix (in the ZYZ convention) 
that takes a free 3- vector of primed components, a', to the frame of reference fixed 
in the laboratory, i.e., a = E{<j), 0,ip) a', its explicit expression is the following [44] : 

/ cos <f> — sin tf> \ / — cos 9 sin 6 \ / cos if — sin i/i \ 
E((f),9,tp) = I sin0 cos0 ( 1 I sin t/j cos ij) I (4j 

V 1 / \ — sin 8 — cos 8 J \ Q 1/ 

» v "• v " v ' 

0(0) <F(V>) 



2 The practice adopted in this section of using different groups of indices, as well as 
different symbols, for the individual coordinates in each one of the sets r, q, e, etc., allows us 
to grasp at a first glimpse the range of values in which each of the indices vary. See table[l]for 
a summary of the indices used, the symbols denoting the sets of coordinates, their meaning 
and the spaces parameterized by them. 
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The position, r' a , of any given atom a in the axes fixed in the system is by construc- 
tion a function, R' Q (w), of only the internal coordinates w, and the transformation 



from the Euclidean coordinates r to the curvilinear coordinates q in (la I may be 
written as follows: 

r Q =R Q (?) =o + E(<I>,8,tP)KW ■ ( 5 ) 

An additional definition of subsets of the coordinates q is motivated by the im- 
position of constraints. Assume that we impose L = N — K independent, holonomic 
and scleronomous constraints to the system: 

a I (q) = 0, I = K + 1,...,N. (6) 

The usual definition [33] is that the constraints are called holonomic if they 
can be written in the form a 1 (q,t) — 0, and scleronomous when the functions 

a T (q) := (a 1 (<?)) i— K+1 do not depend explicitly on time t. They are independent 

if their gradients ((da 1 /dq^)(q)) _ lt I = K + 1, . . . , N, constitute L linearly inde- 
pendent vectors of N components for every point q satisfying the constraints [i.e., for 
every point q such that a 1 (q) = 0, / = K + 1, . . . , N]; or, otherwise stated, the con- 
straints are independent if the matrix with entries (da 1 /dq fl )(q) has range L. In such 
a case, the constraints uniquely define (at least locally) a subspace of W of constant 
dimension N — L = K , called the constrained subspace and denoted by /C. 

Moreover, this independence condition allows, in the vicinity of each point q and 
by virtue of the Implicit Function Theorem [46147] . to (formally) solve eq. ([6| for L 
of the coordinates q, which we arbitrarily place at the end of q, splitting the original 
set as q T = (u T ,d T ), with u T := (u r )^ =l and d T := (d I )f =K+1 . Then, in the vicinity 
of each point q satisfying ([6]), we can express the relations defining the constrained 
subspace, K., parametricallyby 

d z = / 7 (u), I = K + 1,...,N, (7) 

where the functions f T (u) :— (f 1 ( U ))^ K+1 are the ones whose existence the Implicit 
Function Theorem guarantees. The coordinates u are thus termed unconstrained and 
they parameterize JC, whereas the coordinates d are called constrained and their value 
is determined at each point of JC according to Q. As we will see later, in the most 
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Table 1. Symbols and indices used for the different sets of coordinates, as well as the names 
of the spaces parameterized by them. 

Symbol Indices Range Number Name Space 



r a, y8,7, ... 1, ...,n n Atoms 

r,q fj,,v,p,... 1,...,N N — 3n Whole space W 

e A,B,C,... 1,...,6 6 External £ 

w a, b, c, ... 6+l,...,N N — 6 Internal X 

s i,j,k,... 6 + 1, . . . , K M Unconstrained internal E 

d I,J,K,... K+1,...,N L = N-K Constrained internal 

u r, s,t, ... 1,...,K K Unconstrained AC 



general case, the functions / will depend on u, and the constraints will be said to be 
flexible. In the particular case in which all the functions / are constant along /C, the 
constraints are called hard, and all the calculations are considerably simplified. 

Of course, even if AC is regular in all of its points, the particular coordinates d that 
can be solved need not to be the same along the whole spac^] Nevertheless, we will 
assume this to be the case throughout this work, as it is normally implicitly done in 
the literature |48I49|50] . and thus we will consider that K. is parameterized by the 
same subset of unconstrained coordinates u in all of its pointsrj 

Although in general the functions cr{q) in ([6]) may involve all the coordinates q, the 
already mentioned property of invariance of the potential energy function [or the po- 
tential energy plus a term related to the determinant of the whole-space mass-metric 



tensor (see sec. 2.4.3 1] under changes of the external coordinates, e, together with the 
fact (which we shall discuss later) that the potential energy (or the aforementioned 
object) can be regarded as 'producing' the constraints, make physically frequent the 
definition of constraints affecting only the internal coordinates, w: 

a\w) = . (8) 

In such a situation, which will be the one treated in this work, the constrained 
coordinates, d, are contained into the internal ones, w, so that the latter can be split 
as w T = (s T ,d T ), where the first M:=K— 6 = N — L — 6 ones, s T := (s l )f =6+1 , 
are called constrained internal coordinates and parameterize the constrained internal 
subspace, denoted by S. Hence, the constrained subspace can be split as K, = £ x S, 
and the analogue to ^ can be written for this case as 

d 1 = f(s) . (9) 

Finally, if the constraints in ^ are used, together with the Euclidean position 
of any atom may be parameterized with the set of all unconstrained coordinates, u, 



3 One of the simplest examples of this being the circle in R 2 , which is given by f(x, y) := 
x 2 + y 2 — R 2 = 0, an implicit expression whose gradient is non-zero for all (x, y) G AC. 
However, if we try to solve, say, for y in the whole space K,, we will run into trouble at 
y = 0; if we try to solve for x, we will find it to be impossible at x = 0. I.e., the Implicit 
Function Theorem does guarantee that we can solve for some of the original coordinates at 
each regular point of AC, but sometimes the solved coordinate has to be x and sometimes it 
has to be y. 

4 Note that this qualification is unnecessary in the hard case introduced in sec 



2.4.2 



where the uncoupled nature of the functions a(q) makes it possible to solve for the same 
coordinates at every point. 
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as follows: 

r a = TZ a (u) :=K a (e,sJ(s)) = o + E(<f>, 6, <4>)B! a (s, /(a)) = 

=:o + E(4>,e^)Tl' a { S ) , (10) 

where the name of the transformation functions has been changed from R to 1Z, 
and from R' to TV , in order to emphasize that the dependence on the coordinates is 
different between the two cases. 



2.2 Hamiltonian dynamics and statistical mechanics without constraints 

Let us now briefly introduce the formalism needed to tackle the dynamics and sta- 
tistical mechanics of any molecular system such as the one described in the previous 
section. We do this only to fix the notation and to write some expressions that will 
be used later; the interested reader might want to check any of the classical texts in 
the subject, such as [51144152] , 

The central object that determines the behaviour of a classical system is the 
Hamiltonian junction. If one starts by writing the Lagrangian function in terms of 
the Euclidean coordinates r and velocities f (with the over-dot denoting the time 
derivative), 

Wr,r) := \m^) 2 -U{r) , (11) 

where U (r) is the potential energy of the system and the convention prescribing 
summation on repeated indices has been used, then the Euclidean Hamiltonian can 
be easily obtained through the usual Legendre transform: 



1 



ff Eu c(r,7r) := _(7r^ + Z7(r) , (12) 



where 7r are the Euclidean momenta defined as 



dL Euc 



M =l,...,iV, (13) 



with the parentheses around the indices indicating that the convention that prescribes 
summation when indices are repeated is not to be followed. 



The dynamics of the system, once the Hamiltonian in eq. ( 12 ) is known, is given by 
the solutions of a set of coupled, first-order differential equations known as Hamilton's 
equations: 

*, = • <"»> 

If we now want to describe this dynamics in terms of the curvilinear coordinates q 
introduced in sec. |2.1[ which are better adapted to the covalent connectivity of molec- 
ular systems, we can start by taking the time derivative of the change of coordinates 



in eq. ( la) 



(15) 
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and also noticing that, for any given potential energy U{r) expressed as a function of 
r, we can perform the same change of coordinates and define the potential energy as 
a function of q by 

V(q) := U(R(q)) , (16) 

which, as we mentioned, typically only depends on the internal ones, i.e., V(q) = 
V(w). 



If we take these last two expressions to eq. (Ill, we arrive to the Lagrangian in 
curvilinear coordinates, or unconstrained Lagrangian, as we will call it in the rest of 
the manuscript: 

L(q,q):=±q v G vp (q)q>>-V(w) , (17) 

where the whole-space mass-metric tensor (MMT) G vp {q) is defined as 

dRHti dR^jq) 
G vp {q) := -^^^T- ■ (18) 

Following again the usual process of the Legendre transform, we can derive from 



eq. (17 1 the unconstrained Hamiltonian, 



H{q,p):=\p v G v "{q)p p + V[w) , (19) 



where G vp {q) is the inverse of the whole-space MMT in eq. (18) and the canonical 
conjugate momenta are defined as 

p»--=dqT = G ^(<?)<f ■ (20) 

Finally, the Hamilton equations in terms of the curvilinear coordinates and mo- 



menta are formally identical to those in eqs. (14): 



r = , (21a) 

dp,, 

9H(q,p) 

p " = — W~ ' (21b) 

Under the usual assumptions of ergodicity and equal a priori probabilities [52], 
it can be shown that the partition function characterizing the statistical mechanics 
equilibrium of such a dynamics in the canonical ensemble (i.e., at constant volume 
V, number of particles n, and temperature T) is given by 

Z = ^m J e -pH iq , P ) dqdPi (22 ) 

where h is Planck's constant, we denote j3 := 1/RT (per mole energy units are used 
throughout the article, so RT is preferred over kgT) and o:qm is a combinatorial 
number that accounts for quantum indistinguishability and that must be specified 
in each particular case (e.g., for a gas of n indistinguishable particles, ckqm = 
With J dqdp, we denote integration over all positions q^ and their respective momenta 
Pn, for fi = 1, . . . ,N, being the range of integration usually from — oo to oo for the 
momenta, and the appropriate one for each position. For example, if a typical scheme 
for defining the internal coordinates, such as the SASMIC one [JS], is used, bond 
lengths must be integrated from to oo, bond angles from to 7T, and dihedral angles 
from — 7r to 7r. 
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The corresponding equilibrium probability density function (PDF) is thus given by 

P{q > P) = /e ^(^)oW ' (23) 

being P(q,p)AqAp interpreted as the probability of finding the system with positions 
in (q, q + Aq), and momenta in (p,p + Ap), for sufficiently small Aq, Ap. 
In the same sense, the equilibrium average of any observable 0{q,p) is 

(O) = J 0(q,p)P(q,p)dqdp . (24) 

It is also common in the literature to study observables, or properties, which 
depend only on the positions q and not on the momenta p; the native conformation 
of a protein being a notable example of this [5]. In such a case, we can appeal to the 
well-known formula for the iV-dimensional Gaussian integral [53] . 



J-oo loo d "-VdetAf e ' (25j 

where M is an N x N matrix (that must be positive definite in order for the integral 
to be finite), M^ v denotes the entries of M~ l , and v is a (possibly null) TV-tuple, to 
'integrate out' the momenta in the partition function in eq. (22), yielding: 

Z = X (T) f e -P[VW-T§lndetG( q )] dq ^ (2g) 

where 



If we perform the same integration in the joint PDF in eq. (23), the factor xCO 
cancels out and we arrive to the marginal equilibrium PDF in the space of the positions 

q- 

e -/3[V(w)-T§\ndetG(q)] e '0F{q) 
P(?) = J e -p[V( W >)-T§ lndctG(,')] dg , /e-^(9')dg' ' (28) 

where the normal abuse of notation in probability theory has been committed, using 



the same symbol, P, for the two different functions in eqs. (23) and (28), and where 
we have defined 

F(q) := V(w) - TS k (q) , (29a) 
S k (q) := ^lndctG(g) . (29b) 

The notation in the equations above is intentional, in the sense that F(q) can be 
interpreted as a free or effective energy, since it is obtained via the elimination of some 
degrees of freedom (the momenta p) , and it therefore describes the energetics (at least 
the equilibrium one) of the remaining degrees of freedom (the positions q) in a sort 
of mean-field of the ones that have been eliminated. This free or effective character 
is also emphasized by the fact that F(q) depends on the temperature T, even if in a 
simple way, and the analogy can be taken one step further if we regard V(q) as an 
internal energy and the correcting term as a kinetic entropy |54j ; something which is 
compatible with its being linear in RT . 
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Again, the meaning of this last marginal PDF in eq. ( 28 1 is made explicit if we note 
that P(q)Aq is the probability of finding the system with positions in (q, q + Aq), 
for sufficiently small Aq (irrespective of the value of the momenta), and that the 
equilibrium average of any momenta-independent observable 0(q) is given by 



(O) = J 0(q)P(q)dq 



(30) 



Now, as we mentioned in sec. 2.1 the potential energy of molecular systems in 



absence of external fields is typically independent of the external coordinates e, and 
this is why we have written V( w) a nd not V(q) in all the previous expressions. The 
whole-space MMT G(q) in eq. <fT8j) on the other hand, and in particular its deter- 
minant in eqs. (26) and (28), does depend on the external coordinates e. It is thus 



convenient, in the case that we are also dealing with observables 0(w) which are in- 
dependent of the external coordinates, to try to eliminate them from the expressions 
if possible. One can indeed formally always integrate over the external coordinates 
(or any other variables spanning the probability space) in order to get to the corre- 
sponding marginal PDF in the internal space X, i.e., depending only on the internal 
coordinates w. However, until recently, it was not clear if this process could be per- 
formed analytically (specially for the more involved, constrained case in sec. 2.4.1), 
thus yielding manageable final expressions. In a previous work by some of us |55j. we 
settled the issue proving that this is in fact possible and providing the exact analytical 
expressions to be used for the marginal PDF in I. 

To integrate out the external coordinates analytically from det G(q), we only need 
to realize that, since the unconstrained case can be trivially assimilated to a con- 
strained situation with the number of constraints, L — 0, all the results regarding 
the factorization of the induced MMT determinant that we will introduce later in 



2.4.1 apply, and we have the unconstrained analogue of eq. (64) 



det G(q) = sin 2 6» det G'(w) , 



(31) 



with G'(w) the N x N matrix obtained from eq. (65) if the unconstrained internal 
coordinates s are substituted by all the internal coordinates w, and 9 being one of 
the Euler angles introduced in sec. |2.1| 



Now, the external angle 9 can be integrated both in the numerator and in the 
denominator of eq. (28), yielding the marginal PDF in the internal space: 



P(w) 



where we have defined, analogously to eq. (29), 

F(w) := V(w)~TS k (w) 
S k (w) := ^ In det G'H . 



(32) 

(33a) 
(33b) 



Finally, let us remark that every step taken in this section, from the original 
Hamiltonian in eq. (19) to the marginal PDF in eq. (32), is exact. The approximations 
will arrive when we introduce the constrained models in the following sections. 



2.3 Types of constraints 



If we assume that the constraints we shall impose on our physical system can be 
expressed as in Q, where w := (s,d) are typical internal coordinates, such as the 
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ones defined in ref. |45) . then the different types of constraints that can be imposed 
may be classified according to two independent criteria, which, together, constitute a 
definition of the internal constrained subspace S, namely, (a) What coordinates are 
constrained? I.e., what coordinates comprise the set dP. and (b) Do the values of the 
constrained coordinates in S depend on the point s, or not? 



2.3.1 What coordinates are constrained 



Regarding (a) , it is typical in the literature to constrain whole groups of coordinates of 
the same type, classifying them both according to whether they are bond lengths, bond 
angles, etc., and also in terms of the atoms involved in their definition. For example, 
in ref. [8], bond- lengths involving Hydrogen atoms are constrained in a 20 fis MD of 
the Villin Headpiece; in ref. [7j, all bond-lengths are constrained to study the folding 
of a modified version of the same protein; in ref. [55] all bond-lengths as well as bond 
angles of the form H-O-X or H-X-H, being X any atom, are constrained to simulate 
an 80-residue protein; the same scheme of constraints is applied in ref. [57]; and we 
have the constraining of everything but torsion angles in the so-called torsion-angle 
MD |58|59|60|61|62] , 

If the system is close to a local minimum of the potential energy V(w), this 
way of proceeding is based on the intuition that we can associate similar vibrational 
frequencies to internal coordinates of the same type, with little influence of the local 
environment surrounding the atoms that define them. In this way, if we start by 
constraining the highest-frequency coordinates, we remove their vibrations from the 
dynamics, and, as we mentioned in sec. [TJ we can use a larger time-step to integrate 
the equations of motion. This is only intuitive and not rigorous because (i) the idea is 
also used away from minima, where 'vibrational frequencies' are not properly defined, 
(ii) even close to a minimum, the vibrational modes of molecules cannot be directly 
associated to individual internal coordinates, but to linear combinations of them, 
and (iii) even if a vibrational mode is almost purely associated to a given internal 
coordinate, its vibrational frequency does depend on the environment of the atoms 
defining it; the question is how much. 

It is therefore probably a more accurate way of proceeding to define the con- 
strained coordinates in an adaptive manner, depending on the conformation of the 
system; but such schemes either do not exist or they are not very popular in the 
literature as far as we are aware. Therefore, although we remark that any choice of 
constrained coordinates is in principle possible as long as the resulting model is phys- 
ically accurate, we will be thinking here in the simpler case of constraining whole 
groups of internal coordinates of the same type (defining type in the common and 
chemically intuitive way specified above). 

Regarding this point, it is also worth remarking that, although we argued in 
sec. |2.1| that, by virtue of the Implicit Function Theorem, it is locally equivalent to 
express the constraints implicitly, as in eq. (8j, or parametrically, as in eq. the de- 
velopments that can be performed, the final expressions, and the practical algorithms 
derived from them are different in the two cases. We will base all the discussion on 
the formalism in which the coordinates to be constrained, d, can be identified and 
the constraints written as in eq. ([9]), using the unconstrained coordinates, u, to pa- 
rameterize the constrained subspace. In such a case, it is most natural to assume that 
these coordinates, q := (u,d), are the typical ones used in molecular simulation |45j . 
however, the possibility also exists to use the implicit constraints in eq. ([8]) and work 
in Euclidean coordinates. In fact this is the version that is needed for the Lagrange 
multipliers formulation in the whole space and for the implementation of the vast 
majority of practical algorithms, both in the flexible and hard cases discussed later 
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63 49]. Related to this, one can also choose a set of modified curvilinear coordinates 
q in which the constrained ones are given exactly by the a functions in eq. ([8]), and 

therefore the constraints are expressed just as d = d = 0. We will make frequent 
references to these other choices as we advance in the discussion. 



2.3.2 Flexible vs. hard constraints 

Once the constrained coordinates d have been selected, we can ask question (b): Do 
their values in S depend on the point or not? I.e., do they depend on the uncon- 
strained internal coordinates s? We shall call the constraints flexible if the answer is 
'yes' and hard if it is 'no'. 

It is worth mentioning at this point that the wording used to refer to constrained 
systems in the literature is multiple and often misleading. One of the aims of this work 
is to clarify the mathematical definitions behind the words and to provide a consistent 
vocabulary to refer to the different types of models. What we have called flexible 
constraints, for example, are also called flexible in refs. (41(6 4 65 49 50 66] . elastic 
in [57], adiabatic in [68] , and soft in [6914115 70 68 ; whereas the hard case is also 
called hard in refs. [69 41 64 50 68], just constrained in [65] . holonomic in [67168] . rigid 
in [49150] , and fully constrained in [50] . In the light of the more formal discussion in 
what follows, the reader will appreciate that some of these terms are clearly misleading 
{elastic, holonomic or fully constrained), and, in any case, so many names for such 
simple concepts is detrimental for the understanding in the field. Even if the reader 
does not like the choice advocated in this work, she should recognize that a unification 
of the vocabulary is desirable. 

The situation is further complicated by the fact that, when studying the statistical 
mechanics of constrained systems, as we will also see later in great detail, one can think 
about two different models for calculating the equilibrium PDF. The names of these 
models often collide with the ones used for defining the type of constraints applied. On 
the one hand, one can implement the constraints by the use of very steep potentials 
around the constrained subspace; a model often called flexible [71I72I73I74I75I76I77] 
(the worst option due to vocabulary clashes), sometimes called soft [5178] . and some- 
times stiff [79 80 5 48 81 8 2(83] (as we advocate in this work). On the other hand, one 
can assume D'Alembert's principle [44184] and hypothesize that the forces are just the 
ones needed for the system to never leave the constrained subspace during its dynam- 
ical evolution; a model normally called ri gid |71j79|80)72|48)73|74(75|76|85|77|82|83] . 
as we do here, but sometimes also hard [5178165] . which is a bad choice in terms of 
clashes. This interference between the two families of wording is entirely unjustified, 
since, as we will see later, the two types of constraints and the two types of statistical 
mechanics models can be independently combined; one can have either the stiff or 
the rigid model, with cither flexible or hard constraints. 

As we briefly mentioned in sec. |2.1[ the most general case from the mathematical 
point of view is indeed that of flexible constraints, where the functions / in ^ are 
not constant numbers, and it is also the case that arises more naturally from physical 
considerations. 

Indeed, one way of justifying a given constrained model is by comparing it to 
the original, unconstrained one. This choice of the reference is not the only possi- 
ble one; it is also legitimate to compare the constrained models introduced in this 
work directly to experimental results, or to more fundamental theoretical descrip- 
tions such as those based in quantum mechanics. However, we (and many others 
[76I58I85I64I67I59I41I86I87I71I69I88I89I90I77I69I91T5U] ) have preferred to use the clas- 
sical unconstrained model as a reference to be able to assess the influence of the 
imposition of constraints in an incremental way consisting of a series of controlled 
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steps. Although some works have advocated the point of view that the proper jus- 
tification of constrained models should come from a quantum mechanical treatment 
49 54 72 92 56 43 79 83 93 9 J , and the comparison to experiment [H] or to quantum 
mechanics is indeed more relevant in terms of the absolute accuracy of the con- 
strained models, the approximations connecting the former to the latter are many, 
thus obscuring the influence of each individual step. In particular, it is worth men- 
tioning refs. [72 95 73 79 for a clear exposition of the many arguments that ap- 
pear when the comparison is of this second, more convoluted class. This is much 
so that we can find works in the literature that, based on quantum mechanical ar- 



guments, favour the classical stiff model (introduced in sec. 2.4.3) against the rigid 



one (sec. 2.4.1) [72 , or vice versa [35]. Before these discrepancies can be solved and 
a reference classical model can be rigorously derived from quantum mechanics (see 
refs. |96|97|98|99|100|101|102|103|73|81j and also the review by Alvarez-Estrada and 
Calvo in this special issue |104] for the most advanced attempts to do so), we will 
consider the unconstrained dynamics as a more or less 'safe harbor', and as the ref- 
erence to which assess the classical constrained models discussed in this work. It is 
also worth remarking that, in doing so, we are also sidestepping the question about 
the accuracy of the potential energy itself (since we compare to the unconstrained 
dynamics given by the same potential). Although this is an important issue affecting 
the absolute accuracy of the results, and it is very likely that classical force fields are 
not enough to describe some phenomena, e.g., those involving charge transfer [15] . 
we can consider the two sources of error separately and analyze each one of them at 
a time. 

In this spirit, we may think that the fact that allows us to treat some coordinates 
as constrained and some others as unconstrained is that the former need more energy 
to be separated from their minimum-energy values than the latter, and thus this 
separation is more unlikely to occur. In principle, these minimum-energy values as 
well as the energy cost needed to separate the system from them can be defined 
both including the momenta or excluding them, and this can be done either at the 
dynamical or at the statistical mechanics level. 

In this work, we are mainly concerned with the equilibrium properties of con- 
strained models and we point the reader interested in dynamical considerations to 
[fi9l41ll05l67l9ril40ll0r;i73l71l77l82ll07j . as well to refs. [28110811091110] in this same 
special issue. In brief, we want to remark that the approach followed in this work 
is in general simpler and it is subject to less uncertainties than dynamical analyses, 
but, on the other hand, its application must be circumscribed to the computation of 
equilibrium averages. 

In equilibrium statistical mechanics, all the information is contained in the dif- 



ferent partition functions [such as (22)] and PDFs [such as (23)], and all observable 
quantities are calculated through integrals, such as (24). Therefore, it is most natural 
to look for the regions in which the quantity under the integral sign is non- negligible, 
and define our constrained subspace exactly as those regions. Assuming that the 
observables 0(q,p), O(q) are smooth enough functions of their arguments, the con- 
strained subspace so defined corresponds to the values of the constrained coordinates 
that maximize the corresponding PDFs for every value of the rest of variables. Once 
this information is used to define the constrained subspace [i.e., to answer question 
(b)], the issue about how to sample the conformations of the system to produce the 
appropriate equilibrium PDF becomes an algorithmic one (we can use Monte Carlo 
techniques, MD, etc.). If the algorithm used to do the sampling is, as it is often the 
case, a dynamical one, then nothing in the discussed approach guarantees that the 
trajectories obtained are similar in any way to those of the unconstrained system. 
The only thing that we can be sure of is that the equilibrium averages will be in- 
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deed consistent with the approximations involved in the definition of the constrained 
subspace as that of maximum probability. 

In this framework, we have several options to define the most probable values 
of the constrained coordinates d, i.e., the functions / in eq. ^ or (J9J: The most 
immediate one is to look at the equilibrium PDF of the unconstrained system in 



eq. ( 23 ) , before having eliminated the momenta from the description. The values that 
maximize this PDF are those that minimize the unconstrained Hamiltonian for each 
value of the rest of the coordinates, u, and all the momenta, p. This can be formally 
accounted for by stating that the functions / take the values that make H(u,d,p) a 
local minimum for each fixed u and p, i.e., 

H(u,f,p)<H(u,d,p) ,Vd€V(f) , (34) 

where T>(f) is a suitable open set in E L containing the point /, and therefore the 
functions / must satisfy 

f)J-F 

w (u,f,p)=0, I l\ • 1 V . (35) 

with the corresponding Hessian, 

() ~" (uJ, P ), (36) 



dd!dd J 



being a positive definite L x L matrix. 



Now, since the left-hand side of eq. (35) is precisely minus the right-hand side 



of the Hamilton equation in eq. (21b I for the time derivative pi of the canonical 
momentum associated to d 1 , we have that this kind of constraint is equivalent to 
asking that 

Pi = pi(t) =p/(t ) , Vi. (37) 

Using the definition of the momenta in eq. ( p0| , this condition can be seen to 
define the following non-holonomic constraint involving the velocities: 

Gi^qW = Cj , (38) 
being Cj a constant number dependent on the initial conditions. In fact, if we look 



at the implicit definition of the functions / in eq. ( 34 1 or ( 35 1 , we see that they must 
depend on the momenta, i.e., they are f(u,p) and not /(it), which already clashes 
with the holonomic scheme introduced in sec. 12.11 

Although some interpretation problems are associated to non-holonomic con- 
straints [111] , and they are known to not accept in general a closed Hamiltonian 
formalism |112lll3lll4j . they are perfectly valid in the algorithmic sense discussed 
before (i.e., as a tool to define the integration region in which the integrated quan- 
tity is non-negligible in equilibrium statistical mechanics). In fact, when their form is 
linear in the velocities, as it is the case of eq. ( [38] ) , the corresponding equations of mo- 
tion obtained from D'Alembert's principle are very simple |112lll3lll4j . and, even if 
they cannot be expressed as the Hamilton equations of a given Hamiltonian function, 
it has been shown that their equilibrium distribution can be nevertheless computed 
|115j . In [35], for example, this type of flexible constraints are used to build a practical 
algorithm that is shown to be time-reversible and to present good energy conservation 
properties, whereas it is not symplectic; as expected from the non-Hamiltonian char- 
acter of the dynamics. However, it remains to be proved, maybe using the techniques 
in [115) . what is the equilibrium PDF produced by this dynamics. Since in this work 
we are interested in the statistical mechanics of constrained models, we prefer not 
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to deal with this case related to the minimization of H(q,p) until these fundamental 
problems are solved. 

A different option to define the constrained subspace appears if we consider the 
frequent case of momenta-independent observables O(q). In this situation, as we dis- 
cussed before, the momenta can be 'integrated out' from the description, and equi- 



librium averages are computed as in eq. ( 30 1 , using the marginal PDF defined in 



eq. (28). Following the same ideas that we applied in the momenta-dependent case, 



we see that the region of position space in which the integrated object in eq. (30) 
will be non-negligible is the one in which the equilibrium PDF is maximal; or, equiv- 
alently, that in which the quantity F(u,d) := V(s,d) — TS k (u,d) defined in eq. (29) 



is made a local minimum for each value of the unconstrained coordinates u, i.e., we 
ask the functions / that define the constrained subspace to satisfy 

F(u, /(«)) < F(u, d) ,VdE V(f(u)) , (39) 

where 2?(/(it)) is again a suitable open set in R L containing the point f(u), and 
therefore we have that 



ad* [uj{u)) : = i/( s ' /(u) )- T i^^ /! " n = " • ' = K + l v • ,4,1) 

with the corresponding Hessian, 



d 2 F 

dd^dd 1 



(«,/(«)) , (41) 



being a positive definite L x L matrix. 

In contrast with the previous definition, this one (which has never been used in 
the literature as far as we are aware) does produce constraints that are holonomic, 
i.e., the functions / that are now implicitly defined by the equations above depend 
only on position-like degrees of freedom (the unconstrained coordinates u), and the 
constrained subspace is now defined by the holonomic relations in eq. (|7|). Assuming 
D'Alembert's principle, these constraints produce a dynamics which is Hamiltonian, 
which is therefore amenable to symplectic integration methods, and whose equilibrium 
PDF can be easily computed using the same standard ideas that we develop in sec. |2.4| 

In this work, however, for the sake of simplicity, we choose to define the flexible 
constrained subspace in a third, different way which is nevertheless related to both 
of the definitions just discussed. Simply stated, the function of the positions whose 
minimum we use to implicitly define the value / of the constrained coordinates d is 
neither the Hamiltonian H(q,p) nor the free energy F(q), but the potential energy 
V{w), i.e., we ask that, for each fixed value s of the unconstrained coordinates, 

V{s,f(s))<V(s,d) , VdG ©(/(«)) , (42) 

where 2?(/(s)) is a suitable open set in K L containing the point f(s), and we have 
that 

dV 

(*,/(*)) =0, I = K + 1,...,N, (43) 
with the corresponding Hessian, 

d 2 V 



being a positive definite L x L matrix. 

Some points about this last definition are worth remarking before we move on: 
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Fig. 2. Schematic depiction of an unphysical flip due to the use of the global minimum 
criterium for defining the constrained subspace. Three sections of the potential energy surface 
are shown at three different values of x. At xa, the minimum 1 is the global one, whereas, 
at xc, it is 2. At xb, the potential energy of 1 and 2 are the same, and the global minimum 
criterium becomes ambiguous. 



This choice, which is also used in refs. |54|50|79 69 67 88 89 73|41l82)83j and in 
all of the previous works by our group in this topic |48I55I116"| . can immediately 
be seen to come from either of the two previous definitions if some terms are 
neglected. In the momenta-dependent case, we need to neglect the derivative of 
the kinetic energy with respect to the constrained coordinates d in eq. ( 35 1 to 
arrive to eq. ( |43[ ); in the case in which the constrai nts a re defined by minimizing 
the free energy F(q), the derivative of S k (q) in eq. (40 1 needs to be neglected to 



obtain the condition in eq. ( 43 ) . Using that, for any N xN matrix A(x) dependent 
[531 



on a parameter x 



d&etA , „ / x dA 

— = det A ■ tr A — 

Ox \ ox 



(45) 



and also that, trivially, 



A' 1 A 



tr (A- 1 A) = N 



fdA- 1 A 



tr 



A dx ) 



(46) 

the reader can check that the mentioned condition on the kinetic energy implies 
the one on S w {q) (as expected), but the reverse does not hold in general. 
As we will show in the following sections, if suitable internal coordinates are used 
[25], the expression of det G(q) as a function of the coordinates is actually very 
simple. Therefore, not only all the basic ideas in this work are applicable both 
to the definition of constraints using F(q) and the one using V(w), but it is also 
rather straightforward, technically speaking, to add the term associated to the 
minimization of S (q) to the mix. We plan to do this in future works. 
Since, as we mentioned, the potential energy function of molecular systems in 
the absence of external fields does not depend on the external coordinates, the 
constrained values / defined in this last case do not depend on the whole set of 
unconstrained coordinates u, but only on the internal ones s; and the relevant 
constrained subspace is not JC but £ (see sec. 2.1 ). 

For this last case, but also for the previous two, there are at least two reasons why 
the definition based on a local minima is more convenient that a possible definition 
based on a global one (such as the one we erroneously advocated in a previous work 
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by some of us s48J): First, the problem of global minimization is a difficult one; for 
general functions it is known to be NP complete |117I118| . and therefore, even if 
we wished for some reason to define the constrained subspace as the place where 
some function attains its global minimum, there is no guarantee that we will find 
the practical means to compute it. Secondly, such a definition would very likely 
cause unphysical effects: Imagine for the sake of simplicity that we have a system 
with two position-like degrees of freedom, denoted by (x,y), and their respective 
momenta. If we identify the coordinate y as 'stiff' and want to constraint it to its 
global minimum-energy value for each value of the unconstrained coordinate x, 
we could face the situation schematically depicted in fig. [2] In such a case, as we 
smoothly progress along the coordinate x, we can see that the global minimum 
definition might become ambiguous at some value x — xb and then the system can 
instantaneously change the value of its coordinate y to one which is close in energy 
to the previous one but far away in conformational space. The flip is unphysical 
simply because, if the barrier separating the minima is high enough (something we 
shall need in order to consider y stiff), then it will be a very unlikely event that an 
unconstrained trajectory actually performs the flip, and even if it does, it will not 
do it instantly (and with a discontinuity in the constrained coordinates) as in the 
flawed constrained approximation based on global minimization. This situation is 
not academic but actually very common when we have different molecular isomers, 
such as the cis and trans forms of the peptide bond in proteins [BJ. 

Using the definition based on the minimization of V(w) from now on, we can now 



relate the possible answers to question (a) in sec. 2.3.1 to the ideas introduced in this 
section: Following the reasoning that has led us to this point and which is behind the 
definition of the constrained subspace, a selected set of constrained coordinates, d, 
will be a priori a good choice [i.e., a good answer to (a)] if, for 'small' variations Ad 
in d, it happens that V(s, f(s) + Ad) — V(s, /(s)) ^> RT. If this is so, the statistical 
weight, which is proportional to exp ( — V(s, d)/RT), will become negligibly small as 
soon as we separate from the constrained subspace £ by any relevant amount. Albeit 
intuitively appealing (so much that we ourselves argued in this line in a previous work 
48J), we will leave this argument at this point remarking that the scary quotes around 
'small' actually hide real definition problems, as anyone can easily see if she considers 
that, to begin with, some internal coordinates have length units and some others are 
angles. We are at the moment pursuing a more mathematically involved and more 
rigorous definition of the properties that a proper 'stiff' constrained coordinate must 
have, but notice that this only affects the quantitative details pertaining what do we 
mean when we write 'small'; the basic, intuitive idea that the constrained coordinates 
are difficult to change due to energetic reasons is still valid. 

Of course, the answer to question (b) discussed in this section, i.e., whether or 
not (or how much) the equilibrium values f(s) of the constrained coordinates depend 
on the point in £ given by s, is also determined (by definition) by the nature of 
the potential energy used to model the system. Using Quantum Mechanics-based 
potentials, it has already been proved in small peptides that, for a given choice of 
constrained coordinates, their equilibrium values do depend on the conformation s 
significantly [11914811201121] . In this work, we will show en passant, as some previous 
works have already done |49|58)122|41|50j . that the functions f(s) also depend on s 
in the case of the much simpler classical force fields |16ll7ll8ll9l20l2f[ typically used 
for MD of proteins and nucleic acids. This is hardly surprising if we think that the 
force fields potential energy function typically has the following form: 

1 N, No 

Vff(w) :=-J2 Ki a (L - l° a f + ^ E K °« ( « - « ) 2 + y F F (^) + V&W > (47) 

a=l a— 1 
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where l a are bond lengths, 9 a are bond angles, and <p a are dihedral angles. The 
integer Ni is the number of bond lengths, Ng the number of bond angles, and the 
quantities Ki a , Kg a , i° and 9^ are constant numbers. The term denoted by Vp F ((j) a ) 
is a commonly included torsional potential that depends only on the dihedral angles 
^ n the other hand, V FF normally comprises long-range interactions such as 
Coulomb or van der Waals, which involve 0(n 2 ) terms related to all possible atom 
pairs, each one of them depending on the difference of the atomic positions r a — Tp. 

It is precisely this last dependence the one that couples all internal coordinates 
with one another, the one that makes force fields non-trivial, and the one that makes 
the minimum-energy values of, say, bond lengths and bond angles, not the constant 
numbers Z° and 9^ in the force field but some s-dependent quantities. 

Although it becomes clear that the most natural way of thinking about the physi- 
cal origin of constraints presented in this section suggests that they must be (at least in 
some degree) flexible, i.e., that the functions f(s) in eq. (|9| must indeed depend on the 
unconstrained internal coordinates s, the most common approach (by far) in the com- 
putational chemistry literature Il23ll24l75ll^5l78l6lll26l62ll^7l76l85l58l63ll28l5l7^l80l59l90ll^9l86l87l71l77ll30ll3lll32 
has been to consider them as constants, or hard; existing a few works, most of them 
recent, dealing with the flexible case 49 135 65 64 41 69 66 50J, including some previ- 
ous works by us 48 55 116 (although it is worth remarking that the basic idea and the 
fact that it comes more naturally from physical considerations had been mentioned 
as early as 1969 by Go and Scheraga [54]). This is partly due to numerical conve- 
nience but it also might be related to the great influence that classical force fields 



have had in the area. Indeed, if we take a look at the generic expression in eq. (47) 
of such an energy function, we see that the quantities, say, l® a and 0°, appearing in 
the harmonic energy terms associated to bond lengths and bond angles, are readily 
available to be elected as the candidate constant numbers to which equate this coor- 
dinates should we want to consider them as 'stiff' and constrain them; even if they 
are not, as mentioned, the actual minimum-energy values of the bond lengths and 
bond angles in any given conformation of the molecule. If, instead of a force field, we 
consider a less explicit energy function, such as the one produced by the ground-state 
Born-Oppenheimer approximation for the electrons in quantum chemistry [231136] . 
then no candidate number appears before our eyes as the constants to be used in the 
constrained model, and the flexible way of proceeding is even more natural. 

A different but related misconception, probably also stemming from the prevalence 
of force fields in the literature about constraints, is the idea that the steeper the 
potential energy function is when the system leaves the constrained subspace, the 
better approximation the hard one becomes [73)721 . This is indeed true if the potential 



has a form such as the one in eq. (47), because, as the 'spring-constants' Ki a and 
Kg a get larger and larger, the associated harmonic terms overcome the long-range 
interactions, which were, as we argued, the ones responsible for the minimum-energy 
values of bond lengths and bond angles not being the constant numbers 1^ and 9^ 
but conformation-dependent quantities. However, this particular form of the potential 
energy function is not the only possible choice, and one could perfectly argue that 
harmonic terms such as the following ones, which contradict the previously stated 
argument, 



- £ K la (l a - t{s)) 2 + -J2 K o a (*« ~ 0° a (s)) 2 (48) 

OC=l OL — 1 



The terms in Vp F ((j> a ) associated to 'soft' dihedrals, such as the Ramachandran angles in 
proteins, are typically trigonometric functions; whereas those associated to 'stiff' dihedrals, 
such as the peptide-bond dihedral angle uj, are typically harmonic. 
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may model the actual physics of the system more accurately in some cases. Indeed, 
using this kind of functional form, we can have infinitely steep potentials around 
the constrained subspace at the same time that the minimum-energy values of the 
constrained coordinates do depend on s. 

In sec. |2.5| more information is provided about the use of particular energy func- 
tions in the literature and its influence on the flexible vs. hard issue, as well as on 
other approximations. We also point the reader to the discussion about the choice of 
coordinates in sees. |2.4.1| and sec. |2.4.3| which is very much related to the issue of 
the form of the potential energy function in the vicinity of the constrained subspace. 
Finally, in sec. [3j we analyze the difference between flexible and hard constraints 
in a simple methanol molecule, and some works are mentioned in which this same 
comparison has been performed in other practical cases. 



2.4 Statistical Mechanics models 

2.4.1 Rigid model with flexible constraints 

The reasoning based in energetic considerations in the previous section justifies to 
substitute the fact that the system is more likely to evolve close to a given region of 
its internal space, namely, the internal constrained subspace S, by the approximation 
that it evolves only in that region. 

When imposing constraints on a physical system to enforce this approximation, it 
is very common (specially in the theoretical physics literature (44151) but also in com- 
putational chemistry works |5I3| ) to assume not only that the constraints are exactly 
fulfilled at any given time, but also that D'Alembert's principle holds, i.e., that the 
force of constraint does no work (it is entirely orthogonal to the constrained subspace 
£). It can be shown that, assuming these hypotheses, the dynamics of the system, 
described now in terms of the unconstrained coordinates u (i.e., the external ones, 
e, plus the unconstrained internal ones, s) is the solution of the Hamilton equations 
associated to the following rigid Hamiltonian [84 : 

H T {u,rj) := ^r9 rs (u) Vs + V s (s) , (49) 

where the restriction Vs(s) of the potential energy V(w) to the constrained subspace 
£ is defined as 

V E (s):=V(s,f(s)) , (50) 

being the functions f(s) the ones that define the constraints in eq. ^ and whose 
origin has been discussed in the previous section. Of course, in the most general case, 
these functions do depend on s and we are therefore working in the flexible setting. 
Inserting these constraints into the kinetic energy in the unconstrained Lagrangian 



in cq. (17), as well as their time derivatives, 



d = ^ , (51) 
we obtain the rigid kinetic energy: 

K r (u,u) := -u r g rs (u)u s , (52) 
where g rs (u) is the induced mass-metric tensor (or pull-back), given by 
9 „W := OSM + ^<*M + *(»)^M + Bmcf^tm , ,53) 
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being G^ v (u) simply the restriction of the whole-space MMT in eq. (18) to the con- 
strained subspace JC — £ x E, i.e., 



Now eqs. (50) and (52) allows us to construct the rigid Lagrangian, 
L T (u,u) := \u r g rs (u)u s - V E (s) , 



(54) 



(55) 



which finally produces the rigid Hamiltonian in eq. ( 49 ) via the usual Legendre trans- 
form, defining the canonical conjugate momenta as 



T] r : = 



,{u)u s . 



(56) 



The constant-energy dynamics produced by the rigid Hamiltonian in eq. ( 49 ) (via 
the corresponding Hamilton equations) can be implemented in a computer using, for 
example, the algorithms described in refs. |50|49|65] . If constant-temperature dynam- 
ics needs to be performed in this setting, one can use any of the methods available 
for hard rigid simulations [130 63 78 137 , probably after some small adaptation to 
the flexible case |49|138j . See also refs. 109 110|139j in this special issue for more in- 
formation about constrained simulations in ensembles other than the microcanonical 
one. 

With regard to the associated equilibrium statistical mechanics, and making the 
same considerations we made for the unconstrained case in sec. |2.2[ we can write the 
corresponding rigid partition function in the canonical ensemble as 



Z r — 



h K 



-fSH r {u, v ) 



dwdr; 



Hence, the equilibrium PDF is given by 



j e -PH r (u',r,')du'dr)' 



(57) 



(58) 



which has the analogous meaning as the one attributed to eq. (23) in sec. 2.2 being 
the equilibrium average of any observable 0(u,rf): 



{Oh 



0(u,vi)P-[{u,r])&u&ri . 



(59) 



Again, if we are interested in observables that are dependent only on the positions 
u, we can use eq. (25) to integrate out the momenta in eq. (57), yielding 



Xr(T) 



where 



Xr(T) := 



e -p[V S (s)-T§ l»detg(«)] du 



J) 



h K 



(60) 



(61) 



If we perform the same integration in the joint PDF in eq. (58), the factor x^(T) 
cancels out and we arrive to the marginal equilibrium PDF in the space of the positions 



Pt(u) 



e -0[Vs(s)-T§ lndet 9 (u)] 
J e -P[V s (s>)-T§kidetg(u')] du , ' 



(62) 
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being the equilibrium average of any momenta- independent observable 0(u) 



(0) r = / 0(u)P r (u)du 



(63) 



The same invariance of the potential energy of molecular systems with respect 
to the external coordinates, e, that we used in sec. |2.2| to integrate these degrees of 
freedom out is here applicable as well, and this is why we have written Vs(s) and 
not V/c(u) in all the previous expressions. The induced MMT g(u) in eq. (53) on the 



the external coordinates. As in sec. 2.2 



other hand, and in particular its determinant in eqs. (60) and (62), does depend on 

formally integrate over the 



one can always 

external coordinates in order to get to the corresponding marginal PDF in S, i.e., 
depending only on the unconstrained internal coordinates s. However, until recent, 
it was not clear if this process could be performed analytically for general flexible 
constraints, thus yielding manageable final expressions. In a previous work |55j . we 
settled the issue proving that this is in fact possible and providing the exact analytical 
expressions to be used for the marginal PDF in E (this case with constraints is much 
more involved that the unconstrained one mentioned in sec. |2.2[ which had been 
proved long ago by Go and Scheraga using a different development |72)V 

We showed that the determinant of the reduced MMT g(u) can be written as 
follows: 

det g{u) = sin 2 6 dct g'(s) , (64) 

where 9 is one of the external Euler angles and the externals-independent matrix g'(s) 
is given by 



1(3) 






v T (n' a ) 




dH! n . 






dn' a T dn' a 



\ 



(65) 



/ 



where I^ 3 ^ is the 3x3 identity matrix, and the matrices v and J are defined as 



v{n' a ) :-- 



o -z' a y a 

z' a -X' a 

-y a K o 



(66) 



and 




-Ky a -Kz' a 
k 2 +zl 2 -y a z' a 

-y a z' a x'* + y* 



(67) 



respectively, being X' a , y' a and Z' a the three Euclidean components of lZ' a (see sec. 
for a precise definition of these quantities). 



2.1 



This result allows to factorize the exponential in eqs. ( 60 ) and ( 62 ) , and to perform 
a second integration to get rid of the uninteresting external coordinates e, thus arriv- 
ing to the marginal equilibrium PDF in the space parameterized by the unconstrained 
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2.3 



internal coordinates s. i.e., S: 



e -PF T {s) 

Pr(s) = j^ww ' (68) 



where we have defined 



F r { S ):=V s { S )-TS*{s) , (69a) 
S*(s):= -kidet g'(s) . (69b) 

Again, the notation in the equations above is intentional, and we regard F r (s) as 
a free or effective energy, and S*(s) as a sort of a kinetic entropy. 

Using P T (s) above, the equilibrium average of any observable O(s) that depends 
only on s can be calculated as 

(0> r = Jo( S )P r (s)d s . (70) 

Finally, let us remark that the derivation in this section is equivalent to the ones in 
refs. |63)31|69)88|89)79|128|82)83|131| . being the only difference that, in those works, 
different generalized coordinates q are chosen in such a way that the constraint func- 
tions in eq. 

o- 1 (w) := d 1 - f 1 (s) , I = K+1,...,N, (71) 

are among them. Therefore, in these new coordinates q, the constraints and their time 
derivatives are not expressed as in eqs. ^ and (51 ), but simply as 



<F = 0, I = K + 1,...,N , (72a) 
<F = 0, I = K+1,...,N . (72b) 

Despite the earned simplifications that this choice brings, we have preferred to use 
coordinates which do not include the functions a among them because typical inter- 
nal coordinates used in molecular simulations in fact do not include these functions 
[45j . and therefore we can more easily relate to the chemically intuitive quantities 
represented by them (bond lengths, bond angles, etc.) in the formalism. However, 
since many works use the modified coordinates q, we will make a stop now and then 
to mention how this choice affects the final expressions and conclusions in this work. 

The transformation to the new coordinates is actually very simple: 

u r = U r (q) := u r , r = l,...,K, (73a) 
d 1 = D 1 (q) := d 1 - f \s) , I = K+1,...,N, (73b) 

and so it is its inverse: 

u r = U r {q) := u r , r = 1, . . . , K , (74a) 
d 1 = D T (q) := d 1 + f(s) , I = K + 1, . . . , JV . (74b) 

We will see in what follows that the fact that this transformation affects only the 
definition of the constrained coordinates, d, produces remarkable properties and very 
simple rules for most of the transformed objects. 
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This can already be seen in the structure of the Jacobian matrix: 



dQ"(g) 

dq v 



( 1 
Q 



ds i 



(75) 



I J 



and its inverse: 



dq v 



( 1 





1/ 



(76) 



both of which have unit determinant: det J = det J =1. 

It is also worth mentioning that, if the constraints are defined as hard (see the 
next section) , then the a functions in eq. ( 71 ) become 



a 1 (w) := d 1 - 4 , I = K + 1, 



(77) 



being constant numbers, and in such a case the simplicity of the change of coordi- 
nates is even higher, being the Jacobian and its inverse in the expressions above both 
exactly equal to the identity matrix. 



In any case, either if the generalized velocities are related as in eq. ( 51 ) or if they 



are zero as in eq. ( |72b[ ) and in refs. |63|31|69 88 89 79 128 82 83|131j . it is convenient 
to bear in mind that momenta and velocities are only proportional when the matrix 
in the kinetic energy is diagonal (e.g., in Euclidean coordinates), and the statement 
sometimes found in the literature |70|85|73T7 6 128 133 134 about the cancelation of 
the momenta associated to the constrained coordinates is wrong in general. 

Regarding the calculations in this section, the use of the modified coordinates q 
leads to a simpler expression for the induced MMT in eq. ( 53 1 : 



g rs (u) = Gf s (u) , 



(78) 



i.e., it is simply the unconstrained-unconstrain ed b lock of the whole-space MMT G 
evaluated in the constrained subspace (see sec. 2.5 for the implications of this in the 
calculation of Fixman's potential). 
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Moreover, if we use eq. ( 109 ) for the transformation of the whole-space MMT G 
together with the change of coordinates in eq. (74), we can prove that 

g rs (u) = <?£(«) :=G rs (u,0) 



G rs (Q(u,0)) + ^^G Is (Q(u,0)) 



G rJ (Q(u,0)) 



df J (u) 
du s 



G r 



du r 



df'ju) 
du r 



du s 

G Is (uJ(u))+G rJ (u,f(u)) 
df J (u) 



df J (u) 
du s 



du s 



g rs {u) 



(79) 



i.e., the induced MMT is not modified upon the change of coordinates from q to q; it 
is the same as the original one component-wise. 

This directly implies that the calculations associated to the factorization and 
elim ination of the external coordinates in ref. |55j . which led us to expressions ( |64[ ), 
(65), and related quantities, still apply in the new modified coordinates. The only 
change to be done is to place a 'tilde' over the appropriate symbols. 



2.4.2 Rigid model with hard constraints 

As we mentioned in sec. |2.3.2| the most common constrained dynamics considered in 
the literature is not the rigid one with flexible constraints described in the previous 
section, but the one with hard constraints. This option, which is numerically sim- 
pler, and which can be accurate enough for certain applications, is described in the 
following paragraphs. 

The key difference between the two models is that, in the hard case, the con- 
strained coordinates d are not equal to some functions f(s) of the unconstrained 
internal coordinates, as expressed in eq. ([9]), but to some constant numbers do: 

d 1 =d I , I = K +1,...,N . (80) 



This makes all partial derivatives in eqs. (51) and (53) zero, and converts the 
induced MMT g(u) simply in the unconstrained-unconstrained sub-block of the re- 
striction to K, of the whole-space MMT, i.e., 

g rs (u) = G^ s {u) :=G rs (e,s,do) , (81) 

where we have used an over-bar to denote hard objects. This fact has important 



implications regarding the Fixman's compensating potential, introduced in sec. 2.5 



Also notice, that, contrarily to eq. (78), which only held in the flexible case if the 
modified coordinates q were used, we now have this property in the original molecular 
coordinates q. 

Another, more subtle consequence of the use of hard constraints affects the calcu- 
lations of the pa rtial derivatives d7l.' a /ds l , which are a key ingredient of the matrix 



g'(s) in eq. (65). If we realize that the functions 7?/ Q (s) are implicitly defined in 
eq. (FLO) as 

n> a (s):=R' a (s,f(s)) , (82) 
the sought derivatives can be computed using 

^w-^(*/w) + !£(../M)f£c.>- m 
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In the flexible case, all terms in the sum are in principle non-zero, and the fact 
that the derivatives df 1 /ds 1 must be typically found numerically requires the use of 
complicated algorithms to obtain this quantities (see sec.[3]and ref. [116] ). In the hard 
case, on the other hand, all these derivatives are zero and eq. ( 83 ) reduces to 



dW on' 

r( s ) = ^H s ' d °) ' 



d 



(84) 



which contains only explicit derivatives which can typically be computed analytically. 
The efficient calculation of the MMT, MMT determinants, and related quantities in 
refs. [87I124I127I128I86I71I77] . where they consider only hard constraints, is implicitly 
based on this fact. 

Apart from these (important) differences, the dynamics of the rigid model with 
hard constraints is derived following the same steps we took in the previous section, 
and the rigid Hamiltonian in eq. ( 49 1 is still applicable to this case as long as we com- 
pute the induced MMT as indicated in eq. (81). This dynamics can be implemented 
in a computer using many different algorithms (see refs. [130163175] . and references 
therein; and also [140] in this special issue), and, again, under the usual assumptions 
of ergodicity and equal a priori probabilities, its statistical mechanics equilibrium 
can be de scribed in a formally identical way to the flexible case, via eqs. (68), (69a) 
and (|69b|. 



Despite this formal identity, however, it has been shown that the flexible and 
hard dynamics produce dif ferent physical results (see refs. |64I49I69I14T 50 and also 
sec. [3]), and indeed in sec. 2.3.2 we built a strong case for considering the latter as 
an approximation of the former. In sec. [3j we will additionally illustrate how this dif- 
ference in the dynamics produce measurable discrepancies at the level of equilibrium 
statistical mechanics in a methanol molecule. 



2.4.3 Stiff model with flexible constraints 



In the previous two sections, we discussed the dynamics as well as the canonical equi- 
librium of the so-called rigid model, in which not only the constraints are assumed 
to hold exactly but also D'Alembert's principle, which hypothesizes that the com- 
ponents of the forces of constraint tangent to the constrained subspace S are zero 
[44184] , 

Nevertheless, if we accept the unconstrained dynamics as the reference against 
which to assess the accuracy of any approximate model (see the discussion supporting 
this choice in sec. 2.3.2), we must realize that, if the system is in contact with a 
thermal bath at temperature T, even if there exist strong forces that drive it to the 
constrained subspace £, there also exist random fluctuations that will inevitably take 
the system away from this region. Moreover, even in the case in which we can assume 
that the potential energy outside £ tends to infinity (a case in which the constraints 
hold exactly by simple energy conservation arguments |142p . it can be proved that 
the force on S resulting from this limit procedure not only consists of the orthogonal 
component needed to keep the system in the constrained subspace, but also of a non- 
zero tangent component that must be added to the expected gradient of Vs(s), i.e., 
D'Alembert's principle, even in the infinitely steep confining potential case, is not 
satisfied in general |84)96|40)106|143l81j . 

These two in general violated assumptions confirm what we already suspected: 
that the rigid dynamics and statistical mechanics equilibrium are just an approxima- 
tion of the real, unconstrained ones (even if we use the more physically sound flexible 
constraints). Therefore, if the rigid dynamics is to be used in practical simulations in 
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order to gain the advantages discussed in sec. [T] it becomes necessary to assess how 
much error these approximations introduce. 

Of course, in order to compare the rigid statistical mechanics equilibrium PDF 
P T (s) in eq. (68) to the unconstrained one in P(w) in eq. (32) (or the free energies in 



their respective exponents), we should arrive to some objects that are defined on the 
same space, being the most natural candidates P r (s) itself and the marginal uncon- 
strained PDF in the space, S, spanned by the unconstrained internal coordinates s. 
This last function is obtained through the integration of the constrained coordinates 
d as usual: 

(85) 



P{s) := J P(w)dd . 



However, in the general case this is only formal, i.e., for a general potential V(w) 
and general curvilinear coordinates q, this integral cannot be calculated analytically. 
One of the available options is to compute it numerically, with the help of methods 
for calculating free energy differences |144I31I145I32I33] , In this work, on the other 
hand, we have chosen to develop a statistical mechanics model, termed stiff, whose 
effective free energy in E can be analytically obtained and compared to the rigid one 
in a more direct and insightful way than the result of numerically performing the 
integral in eq. ( 85 ) , at the same time that it could be considered an approximation 



to the unconstrained equilibrium containing weaker assumptions than the rigid one. 
Contrarily to the rigid model, the stiff one does account for the statistical weights 
of the conformations close to (but out of) 2J, at the same time that it allows for 
the existence of velocities that are orthogonal to the constrained subspace. Also, 
in the derivation of the stiff model, D'Alembert's principle is never invoked. These 
arguments, together with the knowledge of the controlled steps that take us from 
the unconstrained case to the stiff model, suggest that one could expect the stiff 
equilibrium to be closer to the unconstrained one than the rigid equilibrium not only 
conceptually but also in quantitative terms. Although in this article we will assume 
this to be true and the stiff model will be used as the reference against which to 
assess the accuracy of any constrained model, numerical confirmation will be pursued 
in future works. 

The first step in the derivation of the stiff model consists of Taylor expanding the 
potential energy V(s, d) at a fixed point s, in terms of the constrained coordinates d, 
around the point d — /(s): 



V(s,d) = V(s,f(s)) 



+ 



d 2 V 
dd^dJ ' 



dV . 

dd 1 ^ 
./(*)) 



>/(«)) 



(d 1 f(a)) 



(86) 



(d 1 - f(s))(d J - f J (s)) +0((d- f(s)) 3 ) 



or, realizing that to evaluate in (s, /(s)) is to evaluate in the internal constrained 
subspace S, 



V(s,d) = V s (a) + d I V E (s)(d I - f(s)) 

+ - /'(*)) (d J f(s)) + o((d - f( s )f) , 



(87) 



where we have used the notation for the constrained Hessian in eq. ( 44 ) , and we have 



introduced the notation d{Vs{s) for the gradient of V(w) restricted to S. 
The same expansion can be performed on S^(w), yielding 

S k (s,d) = S k E (s) + d I S ] l J ( S ){d I -f I ( s )) 

+ IdhS^d 1 - f(s)) (d J - /'(*)) + 0((d f(s)f) 



(88) 
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and both of them can be introduced into the integral that appear in the definition of 



the marginal PDF in eq. (85 1 



, r -p[v( W )-TS k (wj\i, 
P{s) := / P(w)dd = . J anr . . . ^ gk , M1 . (89) 

Note that expanding the term associated to the determinant of the whole-space 
MMT det G in the marginal PDF in positions space is in principle more accurate 
that the procedure followed in refs. |79I88I89I63I128I5I90I129I82I83I146| or in our own 
work in ref . [H] , where the dependence of det G on the constrained coordinates d is 
simply neglected, evaluating this determinant directly in d = f(s) and thus effectively 



truncating the expansion in eq. (88) at zero order (often implicitly). In fact, this point 
and the inclusion of flexible constraints in the mix are the only issues that make the 
formalism in this work more general than the ones in refs. [15.79 88 89 82 83]. 

Now, if we assu me tha t the constrained coordinates are indeed 'stiff', i.e., that, 
as mentioned in sec. 2.3.2 for 'small' variations Ad in d, we have V(s,f(s) + Ad) - 
V(s, /(s)) 3> RT, then the quantities under the integral signs, which are propor- 
tional to e - v ( s > d )/ RT ] w iH become negligibly small as soon as we separate from the 
constrained subspace E by any relevant amount. Moreover, we also know that, close 



to E, the Taylor expansions in eqs. (87) and (88) can be truncated at low order. 



Therefore, for each fixed value of the unconstrained coordinates s, we can substitute 
the potential energy and the kinetic entropy terms by their respective low-order Tay- 
lor expansions. Also note that this very same argument can be used if we had chosen 
to define E through the minimization of F (q) instead of V(w). 

In principle, there is no reason to truncate at different orders the potential energy 
V(s,d) and the kinetic entropy S* k (s,(i). 

If we truncate both expansions at order zero, we simply have that the exponential 
does not depend on the constrained coordinates, and the integral J ldd cancels out 



in the numerator and the denominator of eq. (32), yielding the (0,0)-stiff PDF in £ 



P^°\s) = .„„. , (90) 



where we have defined 



F(°^( S ):=V s ( S )-TSl(s) , (91a) 
S%(s) := ^ Indet G' s (s) , (91b) 

being G' s (s) the restriction to E of the matrix G'(w). 

Although this statistical mechanics model is well-defined and may be certainly 



considered for some applications (in sec. 2.5 we mention a number of works in which 
it has been implicitly used), two issues lead to try to improve it: On the one hand, the 
truncation of the Taylor expansions at order zero means that we are only looking at 
points exactly on E, and not 'close' to it, which was one of the weaknesses, we argued, 
of the rigid model. On the other hand, the function we are substituting V(s, d) with, 
i.e., Vs(s), does not have the important property of becoming very large when the 
constrained coordinates get away from E (because it does not depend on them!). 
The fact that we could obtain a finite PDF, Ps°'°\s), is due to the cancellation 
of an infinite quantitjj^J J ldd, in the numerator and the denominator. However, if 

6 At least if any bond length is included among the constrained coordinates, which is 
always the case in practical applications. 
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we consider the partition function of the system after performing the substitution, 

Z = J e~^( Vl) ^~ TSi: ^) dsdd, we see that it is infinite, which is certainly a warning 
and probably a problem. 

If we truncate both expansions at the first order, and we assume that we are 
dealing with flexible constraints defined by the functions f(s) taking the values of the 
constrained coordinates that minimize the potential energy at point s, then the first 
derivatives diVs(s) are zero, and the leading order in the expansion of V(s, d) is again 
Vs(s). Therefore, for the (l,l)-stiff model, the exponent in the integral is proportional 
to V s (s) - TS%(s) - Td I S%(s)(d I - f (s)), and, since the sign of the terms d/S£(s) 
can be anyone, there is no reason for the integral on the constrained coordinates to 
be convergent. If we considered the definition of S involving the minimization not 
of V(w) but of F(q), the situation is slightly different but simple, since in this case 
diVs(s) — TdiSs(s) = 0, and the (l,l)-stiff model would be equivalent to the (0,0) 
one. 

It becomes then clear that the second order terms must be included in the model if 
we want to work with finite quantities. Truncating both expansions at second order, 



and using eq. (251, we find the most general expression for the integral over the 
constrained coordinates that will be considered in this work: 



J exp [ - p(y E - TSl + (djV s - Td I S%){d I - f ) 



+ 1 -(d I -f I )A IJ (d J -f J ))]dd (92) 



2tt\ 



det 2 A exp 



^(d z Vs - Td I S k s )A IJ (djV s - TdjS\) 



where we have omitted the dependence on s, and we have defined 

A I3 := Hfj - TdjjSl , (93) 

in order to lighten the notation. We have also approximated the different ranges of 
integration of the constrained coordinates all by the (— oo,oo) range (notice that 
bond lengths range from to oo, bond angles from to tt, etc.). The reason why 
this is expected to be a good approximation is the same that supports the rest of 
the calculations, namely, that the integrated quantity is only non-negligible close to 
the constrained subspace S, where the constrained coordinates are typically far away 
from the true integration limits. 

Despite this formal result, one needs to notice that, for the above integral to 
exist, the matrix A has to be positive definite. If we had chosen the definition of 
the constrained subspace based on the minimization of F(q), this property would be 
satisfied, since A is precisely the Hessian of F(q) at a minimum. In the case treated 
here, however, the positive definiteness of A is not guaranteed. If w e are in a minimum 



of the potential energy, the constrained Hessian Hfj in eq. (|44| is indeed positive 
definite, but the presence of the term —TdjjS% spoils this property, in general, for 
the matrix A. If such a thing happens, or even if it is 'close' to happening (for 
example, if the smallest eigenvalue of A becomes too close to zero), we must consider 
the possibility of moving to the model based on the minimization of F(q). 

Another issue that may cause the eigenvalues of the matrix A to be small and thus 
render the stiff approximation questionable is that the coordinates we have selected 
as constrained may actually be a bad choice. Indeed, if the statistical weights of the 
conformations do not become significantly small as we move away from S at second 
order in the unconstrained free energy F(w) := V(w) — TS k (w), then we have little 
reasons to believe that the situation might change if we include higher orders. Hence, 
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although the (2,2)-stiff model might not be a very good approximation at every point 



s, we should see the positive dcfiniteness of the matrix A in eq. (93) and the size 
of its eigenvalues as a consistency check, and regard those points at which these 
conditions fail as points in which the (2,2)-stiff model is a bad approximation of the 
unconstrained equilibrium, and the constrained coordinates are probably not 'stiff' 
enough. Although wc have not found conformations in which A (or %; see below) 
contains negative eigenvalues in the practical example in sec. [3j this is probably so 
due to the small size of the molecule studied. For larger systems, we have already seen 
some preliminary indications that these matrices might in fact become not positive 
definite depending on the chosen constrained coordinates. 



Now, using eq. ( 92 ) and recalling that the first derivatives of the potential energy 
restricted to S are zero in the chosen model, we have that the (2,2)-stiff PDF in E 
is given by 

P^)(s) = — , (94) 



ith 



F s W(s) ■= V E (a) - TS k (s) TS A (s) - ^-djS k E (s)A" (s)djS k E (s) , (95) 



2 



being A IJ the entries of A 1 and 



S c A (s) := --lndetA(s) , (96) 

where the letter 'c' in the new entropy stands for conformational, in reference to the 
fact that it appears as the result of eliminating some positions, and not momenta, 
being reminiscent of the conformational or configurational entropies appearing in 
quasi- harmonic analysis |147|90|9Tj . The only difference that we would find if we chose 



the definition based on the minimization of F(q) is that the last term in eq. ( 95 ) would 
not appear. As we mentioned, the (2,2)-stiff model is the most general stiff model in 
this work and it has never been used before in the literature as far as we know. 

The calculation of T-Ls{s), which appears in A(s), can be performed in many dif- 
ferent ways depending on the potential energy we are dealing with (the computation 
of the other part of A is described later). The calculation of det G' IJ (s), on the other 
hand, is entirely 'geometric', depending only on the scheme of internal coordinates 
chosen to describe the system (assuming that the functions /(s) have already been 
calculated). The way to c omp ute detG'(w) analytically has been described before, 
using the analogue of eq. ( [65] ) but without constraints. However, if a specific set of 
internal coordinates is chosen, one can arrive to a much simpler and more explicit 
expression. In our work in ref. |55j , we showed that, if the SASMIC internal coordi- 
nates [IS] for general branched molecules are used (a similar result will hold for any 
well-designed scheme of internal coordinates), we have that 



detc'M = TT < TT £ TT ^ a , (97) 




where l a are the bond lengths, and 9 a the bond angles associated to atom a. The 
masses can be eliminated, since they come out of the logarithm as an additive term, 
which represents a constant change of reference in the free energy, and has no effect 
whatsoever in the equilibrium PDF. Go and Scheraga also proved this result long ago 
using different mathematics [72"] . 
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This expression, together with the additivity of the logarithm, not only helps to 
easily calculate S k (w) [and therefore its restriction S%(s)], but also its derivatives 
with respect to any internal coordinate depending on its type: 

9S k , \ 2R 

■^-(111) = —, (98a) 
dS k 

— (w) = Rcot g e a , (98b) 

orW = 0, (98c) 

where <j) a denotes the dihedral angle associated to atom a. 
Using these results, we can rewrite eq. ([95]) as 



F^\s) := Vz(s) - TS k (s) - TS c A (s) + U^(s) , (99) 

with 

UW(s) ■= -^DifisfiA^&Difis)) , (100) 
being D(d I ) the function defined as 

{2/d 1 if d 1 is a bond length 
cotgd' if d 1 is a bond angle . (101) 
if d 1 is a dihedral angle 

The entries, dj jS k (w), of the matrix whose restriction to E appears in the defi- 
nition of A in eq. (93) can also be easily computed, yielding 

Td 2 IJ S\w) = -^S IJ D 2 {d I ) , (102) 

where Su is the Kronecker's delta, and D2(d I ) is the function defined as 

f2/(d 1 ) 2 if d 1 is a bond length 
1/ sin 2 d 1 if d 1 is a bond angle , (103) 
if d 1 is a dihedral angle 



allowing to rewrite eq. ( 93 1 as 

AlJ ( s ) := Hfj(s) + - p 8 IJ D 2 {f I {s)) . (104) 

Note that both D(d I ) and Z?2(d 7 ) are singular only in points that have little 
physical sense and that will be never reached in a practical simulation: l a = or 
6 a = 0, 7r. Not only the energy will be very large (even infinite) at these points, 
but also the change of coordinates in eq. ([!]) from the Euclidean to the curvilinear 
coordinates becomes singular. 

Apart from the (2,2)-stiff model, we can also define two other models which yield 
convergent integrals for any point s. However, they are derived truncating the poten- 
tial energy and the kinetic entropy at different orders; something which is, in principle, 
not justified. 
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Particularizing eq. (92), we have that the (2,0)-stiff effective free energy is given 

by 

F^( S ) := V s (s) - TSl{s) - TS< n (s) , (105) 

with 

Sh(*) -=-j^detU s {s) . (106) 

As we discuss in sec. |2.5| this model has been implicitly used in some works, 
including one by some us |48) . 

The (2,l)-stiff free energy, in turn, includes an additional term: 

F S ^( S ) := V s (s) - TS k s (s) - TS C H ( S ) + U^(s) , (107) 

defined as 

UW(s) := -~D(d I (s))U I s J (s)D(d J (s)) , (108) 

where we have used again that the first derivatives of the potential energy are zero 
in E, and H 1 ^ (s) are the entries of the inverse matrix of the constrained Hessian, as 
usual. If the definition of the constrained subspace based in the minimization of F(q) 
is used, then the (2,0)- and (2,l)-stiff models become identical. 

Now, if the modified coordinates q defined in sec. |2. 4.1] are used in the calculations 
in this section, several points need to be re-examined. First of all, we have the following 
relationship between the whole-space MMT in the coordinates q and the one in the 
coordinates q (indeed, we have it for any 2-times covariant tensor 46J): 

Hence, if we use the fact that the determinant of a product of square matrices is 
the product of the determinants, and also the property mentioned in sec. 2.4.1 about 
the Jacobian determinant of the change of coordinates being 1, we have that 

detG(g) = detG(Q(q)) . (110) 

Note that this not only guarantees the factorization of the external coordinates 
in eq. (31), but it also allows us to use the explicit expression (97) in terms of the 
molecular coordinates q to compute det G{q) as long as we insert the transformation 
functions Q{q) into the appropriate functional places corresponding to bond lengths 
and bond angles. 

Also using the structure of the Jacobian in cq. ( 75 1 , it is easy to prove that the 
transformed potential energy, 

V(q) := V(Q(q)) , (111) 
and indeed any scalar function (such as detG above, and hence S k ), satisfies 

dV dV , 

as well as 

d 2 V d 2 V , 

These transformation rules guarantee that all the final expressions in this section 
are identical in the case that we wished to use the modified coordinates q, as long 
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as we 'place tildes' over all relevant objects, such as Ti(q), A(q) or S k (q). But not 
only that, they also allow us to compute the objects in the original coordinates q 
and then simply perform the substitution Q(q) to find the transformed quantities. 
This is very convenient, since the coordinates q are defined in terms of the functions 
/(s), which are typically impossible to express analytically, while the coordinates q 
are simple functions of the Euclidean coordinates, thus admitting compact algorithms 
for computing quantities such as the Hessian of the potential energy, or the functions 
D(d) and D2{d) in eqs. (101) and (103), respectively. 

It is also worth remarking that, in the derivation of the stiff models in the coordi- 
nates q, the Taylor expansion of the potential energy in eq. ( 86 ) would be performed 
at d — and take the simpler form 



V(s,d) = V(s,0) + 



dV 



d 2 V 



dd J dd J 



d J d J + 0(d 3 ) , (114) 



which allows to make contact with some of the potential energy functions in previous 
works [67 88 89 79 83 69] , realizing that the intention of the authors may not be to 
restrict themselves to a particular case of the completely general potential considered 
here, but just that they are using the modified coordinates q, instead of the molecular 
coordinates q used in this work. See however the discussion in sec. |2.5| about the 
neglection of the dependence of the Hessian on s, which is not a general feature of 
cq. (114), irrespective of the coordinates used. 

Finally, to close this section, let us raise an important and subtle point related 
to the stiff model and the dynamics of the system: As we have emphasized in pre- 
vious sections, the analysis presented in this work is completely at the equilibrium 
statistical mechanics level; all approximations are performed in equilibrium PDFs 
and all constrained models are defined in these terms too. In the previous rigid sec- 
tions [2AT] and [2A2J we briefly commented that the rigid equilibrium can be sampled 
using available dynamical algorithms such as SHAKE [130 63 140 or any of its more 
modern descendants (see ref. [78] and references therein) for the hard case; and the 
algorithms in refs. [50149165] if the constraints are flexible. That this sampling can be 
achieved by such algorithms is easily understandable if we recall that the rigid model 
is obtained from a Hamiltonian [the one in eq. (49)] and thus the corresponding 
discretized Hamilton equations (properly thermostated |109| ) will do the trick. More- 
over, the fact that the fastest vibrations are absent from this dynamics makes the 
integration algorithms potentially more efficient than the unconstrained dynamics, 
and the whole scheme is thus useful. 

In the case of the stiff model introduced in this section, the situation is entirely 
different. We may of course ask the theoretical question of which would be a dynam- 
ics that directly samples the stiff statistical mechanics equilibrium, and we may even 
find some dynamics that do the job. However, any naive version of such a dynamics 
will in principle contain the fast vibrations that we sought to eliminate with the use 
of constraints, and therefore any algorithm that implements it would probably be 
as numerically costly as the whole unconstrained dynamics (at least). Since the only 
rationale behind the development of the stiff model was to provide a good approx- 
imation to the unconstrained system in order to be able to assess the accuracy of 
the rigid equilibrium, using a new 'stiff' dynamics instead of the unconstrained one 
itself seems a rather unprofitable (and convoluted) enterprise. In principle, the stiff 
model should be regarded as an approximation to the unconstrained equilibrium at 
the level of statistical mechanics, not at the level of the dynamics. However, as we 
will discuss in sec. |2.5[ there is a way of sampling the stiff equilibrium using a rigid 
dynamics with an additional correcting term (see fig. [3] for a schematic summary of 
these ideas). Whether or not the calculation of this correcting term can be performed 
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fast enough so that we end with a dynamics which has the accuracy of the stiff model 
while retaining the efficiency of the rigid one is, however, a topic that we will not 
analyze in depth in this work. 



2.4.4 Stiff model with hard constraints 

In principle, we can also consider the stiff model defined using the hard constraints in 
The main difference with the flexible case appears in the Taylor expansions 
which become 




V[a, d) = V s (s) + djVsis^d 1 - 4) (f 15) 

+ \nfM(d x -4)(d J -d{) + o((d-d ) 3 ) , 



and 



5 k ( s , d) = S k s (s) + djSUsXd 1 - 4) (f 16) 

+ l&vS^Wd 1 - 4){d J - 4) + 0((d d ) 3 ) , 

where we have used again an over-bar to distinguish objects associated to the hard 
constraints from their flexible counterparts (meaning they are evaluated at d = do 
instead of d = f(s)). 

Also, since the values d = do do not need to minimize the potential energy 
at the point s [nor the free energy F(q)], the derivatives diVs{s), diS%{s) and 
diVs(s) — TdiS^j(s) are all non-zero in general, and neither the Hessian / H E {s) nor 
the associated matrix A(s) are guaranteed to be positive definite in the hard case. 

All calculations for the hard stiff models are formally identical to the ones in the 
previous section as long as we include the terms diVs(s) and 9j5^(s), and substitute 
the functions /(s) by the constant values do. 

In particular, the hard (2,2) -stiff model has free energy 

F^ 2 Hs) := V s (s) - TS%(s) - TS\{s) + U^(s) , (117) 

with 

S%(s) := | In detG'^s) , (118a) 

St{s) := --lndetl(s) , (118b) 



and 



2 



U^H-s) ■= -\ (djVzis) A"(s) (djV s (s) 



being 

As in the flexible case, the substitution of the second order expansion of the free 



Au(s) := nfj(s) + ^6 IJ D 2 (d I ) . (119) 



energy into the integral in eq. (89) is only justified if A is positive definite and has 
large enough eigenvalues. In the case of the hard (2,2)-stiff model, the failure of such 
a property can be caused not only by the fact that the chosen constrained coordinates 
are actually not very 'stiff' or because we have decided to minimize V{w) instead of 
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F(q) (as it happened for the flexible case in the previous section), but also because of 
having expanded around a value, do, which is not a minimum of the potential energy 
[or of F(q)], thereby possibly distorting the observed curvature encoded in H and 
A. 

The only hard stiff model that presents an equilibrium PDF that is convergent for 
all points s (but a partition function that is infinite) is the (0,0)-stiff model, whose 
free energy is given by 

Fj°'V(s):=V E (s)-TS%( S ) . (120) 

As in the flexible case, the first order truncation produces models that are, in gen- 
eral, non-convergent, while we can still define two models that may be well-behaved 
depending on the properties of the hard constrained Hessian H s : the (2,0)- and (2,1)- 
stiff models. Notice, on the one hand, that these two models were convergent at all 
points s in the flexible case, because the flexible constrained Hessian is, by construc- 
tion, positive definite. On the other hand, the same caution applies in the hard case 
regarding the a priori justification of truncating the potential energy V (s, d) and the 
kinetic entropy S k (s,d) at different orders. 

The hard (2,0)-stiff model has free energy 

Fi 2 -°)( S ) := V s (s) - TS k s (s) - TS c n (s) + U^(s) , (121) 

with 

UM(s) := -\d I V s {s)U I s J {s)d J V s {s) , (122) 

and 

Sui s ) '■= --lndetT^s) , (123) 
whereas the free energy of the hard (2,1 )-stiff model is given by 

Fi 2 '°)( s ) := V s (s) - TS k s (s) ~ TS c n (s) + U^(s) , (124) 

being 

U^(s) := -\ (dfoW ^f-) H¥(s) (djV s (s) - . (125) 

It is also worth noting here that, in the hard case treated in this section, there is 
no difference between minimizing V(w) or F(q), because the constrained subspace is 
not defined through any minimization, but it is postulated to simply be d = do. 

Some of these models have been implicitly used in the literature before, as we 
mention in the next section. In fact, being the hard constraints more popular, their 
use have been more common than that of their flexible counterparts in sec. |2.4.3| 

2.5 Summary, comparisons and the Fixman's potential 

In fig.|3j a schematic depiction of the models discussed in the previous sections can be 
found, together with the relationships among them, as well as the basic mathematical 
expressions used to describe their respective statistical mechanics in the constrained 
subspace. In the figure, only the flexible versions of the constrained models have 
been detailed, and only the more accurate (2,2)-type of stiff model is shown. An 
equivalent scheme will be obtained if we define the constrained subspace S using the 
hard relations d = do or we truncate the Taylor expansions defining the stiff model 
at a different order than 2. 
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Fig. 3. Summary and relationships among the models discussed in this work. For a precise 
definition of the different mathematical objects and the meaning of the models, see the main 
body of the article. 



In words, we have the unconstrained dynamics, which, we have argued, is the 
reference to which we will judge the accuracy of any constrained model. On the other 
end of the scheme, we have the rigid dynamics, which is potentially faster than the 
unconstrained one because the fastest vibrations have been eliminated from the sys- 
tem by constraining the associated degrees of freedom. Despite this better numerical 
profile, the rigid dynamics, and its corresponding statistical mechanics equilibrium 
in the constrained subspace S are based on some questionable assumptions, such 
as D'Alembert's principle, or the imposition that the dynamics takes place exactly 
on S, thus forbidding the system to present velocities that are orthogonal to this 
space. These strong approximations require that we assess the accuracy of the rigid 
model against the reference unconstrained situation. In this work, we choose to do 
so at the level of equilibrium statistical mechanics, where the key object is the PDF 
and the free energy that appears in its exponent. Given the typical complexity of 
molecular potential energy functions, the integral needed to marginalize the uncon- 
strained PDF to the constrained subspace £ and thus be able to compare it with 
the rigid one is too complicated to be performed analytically. This is why we have 
introduced an approximation to the unconstrained equilibrium, which has S as its 
natural probability space, which is based in weaker assumptions than the rigid model, 
and which is termed stiff model. Now, by comparing the rigid equilibrium to the stiff 
one, the accuracy of the former can be assessed. That they are indeed apprecia- 
bly different has long be known in the literature and many works have discussed it 
[75l7HI85l88l89ll28l5l54l72ll48l7HI8()l41l79l86l87l71ll49l6HI77l82l8HliaHI81j . 

The most natural way to perform this assessment is to compare the two effective 
free energies that appear in the exponents of their respective PDFs. For example, if 
we consider the more general case of flexible constraints in sees. 2.4.1 and 2.4.3 and 
we use the most accurate (2,2)-stiff model, we can define the following function of the 
unconstrained coordinates s: 



V F (s) := F^(s) - FAs) = TS^(s) - TS k s (s) - TS c A (s) + U^(s) , (126) 
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which quantifies the difference between the two free energies, and contains many of 
the terms related to the MMT and Hessian determinants that we have introduced 
in the previous sections (but notice that the constrained potential energy Vs(s) has 
been canceled out). 

This function was first introduced by Fixman |128j (although in a much more sim- 
plified version) and it is thus normally called Fixman's potential. Apart from quan- 
tifying the discrepancies between the equilibrium statistical mechanics of the rigid 
model and the, in principle, more accurate stiff one, it can be easily seen that the ad- 
dition of Vf(s) to the original potential energy Vs(s) in the rigid dynamics produces 
an evolution of the system whose equilibrium is now given by the stiff PDF, at the 
same time that the fastest vibrations remain absent from the dynamics (assuming 
that V f(s) has not reintroduced them). Indeed, if we modify the rigid Hamiltonian 



in eq. (49) and use the following one instead: 

H?(u, n) := lvr9 rs (u) Vs + V E (s) + V F (s) , (127) 

it is immediate to see that the canonical equilibrium PDF of the corresponding dynam- 
ics (after integrating out the momenta and the external coordinates as we described 
in the previous sections) will be given by 

e -f>F?(s) 

* {s) = T^TMoV ' (128) 

where the obtained free energy Ff (s) can be checked to be exactly the stiff one after 
some simple cancelations: 

F?{s) = V E (s)-TS*(s) + V F (s) 

= V s (s) TS*(s) + TS*(s) - TS%(s) - TS c A (s) + U^( S ) 
= V s (s)-TS%( S )-TS A (s) + U^(s) 

=:F^(s) . (129) 

This practical use of the Fixman's potential to sample the stiff equilibrium using a 
modified rigid dynamics has often been discussed in the literature [131186114918918817717118517311051761851801411771691155] , 
and it is probable the main motivation behind its first introduction |128j . In this work, 
however, no dynamical treatment is going to be performed, and we will use Vp(s) sim- 
ply as a way of quantifying the difference between the stiff and rigid models at the 
level of equilibrium statistical mechanics. As we discussed before, this comparison can 
be understood as an assessment of the accuracy of the rigid equilibrium and, indi- 
rectly, as an evaluation of how much error the approximations behind the rigid model 
(including D'Alembert's principle) produce in the resulting statistical mechanics of 
the system. 

We also point out that, if the justification of constrained models is not based in 
their comparison to the classical unconstrained case, as in this work, but to quantum 
mechanical calculations, then it can be argued that the equilibrium produced by the 
rigid dynamics is the 'correct' one, and therefore no corrections such as Vp(s) are 



needed [35] . We discussed this issue in detail in sec. 2.3.2 

Apart from the introduction of Vf(s), Fixman proved in his celebrated paper 
128 a theorem that, as we will see later, allows to simplify some MMT determi- 
nants in some cases. This useful result has been used in some works to compute 
certain PDFs even when no comparisons between stiff and rigid models are sought 
[144I63I3T] . Additionally, the issue of the comparison between constrained models in 
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which the orthogonal velocities to the constrained subspace are activated and those in 
which they are assumed to be zero is central to the interpretation of the methods for 
computing free energy differences along a reaction coordinate |144I31I1 50 32 33 146] ; 
however, despite the formal similarity, the reasons for imposing constraints and the 
physical considerations behind them are very different in the two cases. We suspect 
that these facts add even more confusion and complexity to some of the accounts 
about the topic. 

It is also worth pointing out that the dynamics produced by the modified Hamil- 



tonian in eq. (127) needs not to bear any resemblance to either the original rigid 
dynamics or the unconstrained one. The only convenient features of this dynamics, 
by construction, are the ones we have mentioned: (i) its equilibrium statistical me- 
chanics is the stiff one, and (ii) assuming that Vf(s) does not introduce new fast 
oscillations, the dynamics is still rigid-like, in the sense that the fastest degrees of 
freedom have been removed and thus it can be integrated with a larger time-step 
than the unconstrained dynamics. In |40| . for example, it is shown, under reasonable 
assumptions, that the averaged dynamics in the infinitely stiff case satisfies some 
modified equations of motion that evolve the trajectory in the constrained subspace, 
and whose modification is not the force corresponding to the Fixman's potential 
introduced here. This discrepancy could have several origins: 

— The differences in the dynamics could be averaged out at equilibrium and the 
resulting PDF may turn out to be the same for the two cases [as we said, the 



dynamics produced by eq. (127) are only good for statistical mechanics sampling]. 
Since the dynamics in eq. (12) of [40] is, in principle, non-Hamiltonian, the tech- 
niques in [115] should be used to obtain its equilibrium PDF. 
As we mentioned in sec. [T] and we will develop further later in this section, the 



calculations in this work and the definition of Fixman's potential in eq. ( 126 ) 
are the most general in the literature as far as we are aware. In every other 
work, including [40j . the version of the Fixman's potential that is used includes 
approximations. Therefore, a comparison of the modified dynamics in |40| to the 



one produced by the Hamiltonian in eq. (127) could yield different conclusions 



— Although the truncation of the Taylor expansion that took us to our definition of 
the stiff model in sec. |2.4.3| seems very much related to the 'infinitely stiff' limit 
of the equations of motion in [30], there may exist subtle differences that could 
make the two physical situations sit on a different basis. 

We do not pretend to solve this issue here, and the only reason why we have 
commented on it is to highlight the fact that dynamical analyses of constrained models 
are much more involved than the straightforward approach in this work, based in 
equilibrium properties. If one disregards the accuracy of the trajectories, the fact 



that the modified Hamiltonian in eq. (127) produces the stiff equilibrium through a 



constrained dynamics is very easy to prove and uncontroversial, as we have shown. 
Since the stiff model does not invoke in its derivation the D'Alembert principle, it is 
also clear that Vp(s) is an appropriate way of quantifying the accuracy of such an 
approximation; even if it were not the appropriate choice for doing the same thing at 
the dynamical level. It is also worth remarking that, doing similar dynamical analyses 
as those in |40j , Reich finds a correcting term for the unconstrained coordinates which 
does seem to agree with the Fixman potential |141I69] , 

Having these issues in mind, as we mentioned, the calculations in this work are the 
most general in the literature as far as we are aware. This is so because of a number 
of points: 



— Apart from the requirement that they are 'stiff', no special property has been 
assumed about the constrained coordinates d. On the other hand, in many works 



Will be inserted by the editor 



39 



63|69 67 88 8 9|128|79)82|83)131j . the constraint functions a 1 (w) := d 1 - f T (s) 



in eq. (71 1 have been chosen as the set d. In such a case, as we have discussed 
in the previous sections, many simplifications occur. For example, the induced 
MMT g in the rigid models in sees. |2.4.T] and |2"A2| is actually the unconstrained- 
unconstrained sub-block of the whole-space MMT restricted to /C, G/c- Therefore, 
the theorem by Fixman |128| applies, and we have the following relation, which 
helps to more efficiently compute the difference TS^(s) — TS^j(s) in Vp(s), and 
which has been used in many previous works [TIT 76 63 88 89 5 179186] : 

detg'(s) det g(u) A f , , iqn s 
a a 4- i \ = det h K.(u) , (130) 



(131) 



detC^s) detG K {u) 

where h denotes the constrained-constrained sub-block of G _1 : 

hI j ^dQ 1 1 dQ J 
dr 11 dr^ 

The reason for expressing the above quotient of determinants as a function of det h 
is that the entries of h are in principle very easy to calculate if q are the typical 
molecular coordinates consisting of bond lengths, bond angles, and dihedral angles 
[45] , since each one of the functions Q depends only on a few Euclidean coordinates 
r. This is however not true if the modified coordinates q including the a functions 
are used. Indeed, if we compute the matrix h in the q coordinates and use the 



transformation rules in eq. ( 73 ) , we find that 



= g(«»)f£w»)£f£wfl)f£(««) 

_ df 1 dQ r 1 dQ s df J _ dQ 1 1 dQ s df J 
du r dr^ nifi dr^ du s dr^ dr^ du s 
_ df 1 dQ r 1 dQ J + dQ 1 1 dQ J 
du r dr^ dr^ dr^ dr^ 

= hIsd ^ - ^ hrJ + h " . ( 132 ) 

ou r ou s au s ou r 



where we have extended the definition in eq. (1311 to run in the whole range 
jj,, v = 1, . . . , N, and we have dropped the dependencies towards the end of the 
calculation to make the expressions lighter. Now, we can conclude that, even if 
in the flexible case one can use the modified coordinates q and the theorem by 



Fixman in eq. ( 130 ) still applies, it is not very useful. This is so because, although 
the terms related to the untransformed matrix h in the above expression are still 
easy to compute, the derivatives df/du are typically non-trivial, and they must 
be found numerically (see sec. [3] and ref. |116| ). In the hard case, however, we have 
f I (s) = dp, and the relation between h and h becomes 

h IJ (q)=h IJ (Q(q)) , (133) 



I.e., in the hard case, eq. ( |130[ ) is not only true, but also useful. In sec. II of 
ref. [7S] , a very simple model illustrates how the choice of coordinates discussed in 
this point might have implications on the whole formalism apart from the Fixman 
theorem. In sec. 2.4 of ref. [57], the authors briefly agree with us about the fact 
that the use of the coordinates q assumes that they have already been identified 
(computed), but this is non-trivial in general. 
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We have worked in the more general flexible setting, which includes the more re- 
stricted (and much more popular |1^3ll^4l75ll25l78l6lll26l62ll27l76l85l58l63ll28l5l72l80l59l90ll29l86l87l71l77ll30ll3 



hard models as a particular case. Note that, according to the discussion in sec. 2.4.3 
it is clear that this issue is coupled to the choice of coordinates discussed in the 
first point of this list. Indeed, if the modified coordinates q introduced in the pre- 
vious sections are used, we have that the Taylor ex pansi on of the potential energy 
around the constrained subspace is given by eq. (114 1. Since the orders higher 
than 2 will also be multiplied by the undisplaced constrained coordinates d, we 
have that, in the coordinates q, the constraints defined through the minimization 
of the potential energy arc effectively hard, i.e., d = 0. However, despite this re- 
formulation, the vast majority of the works in the literature do not refer to this, 
but actually to constraints with the form d — do being the coordinates d typically 
bond angles and bond lengths, which is not a reformulation but a restriction to a 
more particular case than the general, flexible one discussed here. 
The Taylor expansion used to build the stiff model in sec. |2.4.3| has been truncated 
at the highest possible order compatible with the possibility of having analytical 
integrals, thus producing more accurate corrections, such as the ones including 
the derivatives of the determinant of the whole-space MMT G. This Taylor ex- 
pansion of the term —TS k (q) has not be accounted for in any previous work as 
far as we are aware (even in those works in which the treatment is most general 
[48 79 88 89 82 83 ), thus effectively making the most accurate stiff model consid- 



ered so far the (2,0) one in sees. 2.4.3 and 2.4.4 Also, it is often assumed in the 
literature that the zero-order restriction to S, —TS^(s), does not depend on s 
|75I72I90I129I87I133I134] . The first approximation is wrong in general both in the 
flexible and hard cases [it is sufficient to look at eq. (97), where the dependence 
of det G' on the typically constrained bond lengths ana bond angles is explicit] . 
However, if hard constraints are chosen, and all bond lengths and bond angles 
are constrained, then det G' is a constant and it can be canceled out from the ex- 
pressions [751721129187171177] . On the other hand, if at least some bond angles are 
considered unconstrained (as it is common in the literature [76I8I7I56I57I133] ). we 
have that (i) —TS k (q) depends on the constrained coordinates, thus suggesting 



the necessity of performing the Taylor expansion in eq. ( 88 1 , and (ii) even the 
zero-order restriction to E, —TS ] ^(s), does depend on s, also in the hard case. 
The Taylor expansion mentioned in the previous point has also been performed 
here on the potential energy, thus producing correcting terms which include the 
Hessian matrix. In fact, since the best stiff model that appears in the previ- 
ous literat ure is, as we mentioned, the (2,0) one, the matrix A in sees. |2.4.3 
and 2.4.4 is never found, only the Hessian (at most). This correction related 



to the Hessian has not been considered in a number of works, thus reducing 
the stiff model used to the (0,0)-type |75l85l63ll28l5l80l90ll3lllS3TTM] (some- 
times legitimately so, because of the restricted form of the potential energy used 
|71|86|87|77j ). and it has been included (or at least mentioned) by some others 
17616718818915417214911291791821831691 . However, even in the second family of pa- 
pers, it is often assumed that the Hessian does not depend on the unconstrained 
coordinates s [76 72 129 69 , which leads to effectively discard it. It is also worth 
pointing out that this last approximation is sometimes justified by the fact that 
the force field- like potential energy functions [such as the one in eq. (47)] typ- 
ically have conformation-independent spring constants 76 72]. However, even if 
we accept to particularize the calculations to classical force fields, the spring con- 
stants can only be assimilated to the Hessian in the infinite stiffness limit, which 
is not clear to be reached in force fields, specially for bond angles. Finally, we can 
think that the assumption that the zero-point energy of the quantum harmonic 
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oscillators associated to the constrained coordinates is s-independent in quantum 
models is somehow related to this point [54172149] , 

The potential energy used to model the system and to define the constrained 
subspace in the flexible case has been considered to be completely general. It is 
common in the literature to explicitly or implicitly assume that a force field-like 
function [such as the one in eq. J 47|] is used [69l76l67l88l89ll28l54l72l41l79l83j . 
or even simpler ones [7l 86 87 77J, which leads to more restricted versions of the 
objects computed here. For example, in ref. [71] . (and similarly in [77]) a polymer 



is considered whose potential energy is just the harmonic terms in eq. (47), i.e 



V(w) 



(134) 



In such a case, the flexible constraints, as defined in this work, become d = do, 
i.e., exactly hard; with the associated simplifications discussed in the previous 
points. In the celebrated paper by Fixman [128 a similar choice is made. Also, 
as we mentioned in sec. |2.4.3[ the choice of coordinates discussed in the first 
point of this list is related to the issue of the form potential energy. The Taylor 
expansion of the potential ene rgy a round the constrained subspace in the modified 



coordinates q is given by eq. (114), which is in agreement with the potential used 
in some previous works [67 88 89 54 41 79 83 , implying that maybe the authors 
were not considering a restriction to a force field-like function, but simply using the 
coordinates q. Remember however what we discussed about these coordinates in 
the first point, and also note that in some of these studies it is typically assumed 
that the Hessian is independent from the conformation 41 69j, which is not a 



general feature of eq. (114). In sec. II of ref. [75], the relationship between the 
coordinates chosen and the form of the potential energy is clearly illustrated using 
a simple 2-dimensional system, and see also sec. |2.3.2"] for more implications about 
this point in relation to flexible constraints. 

In summary, in this work, all quantities appearing in the expressions: the MMT 
and Hessian determinants, the functions /, the derivatives of the potential energy, 
etc., have been considered to change with every variable that they can in principle 
depend on. No 'constancy' assumption has been made in our development. They have 
been expanded to the highest possible order compatible with analytical results, and 
no special form of the potential energy has been assumed. 

The first issue in the above list affects the form of the final expressions and there- 
fore complicates the 'translation' between our results and those in other works, but 
it does not introduce additional approximations. The rest of the assumptions (that 
we did not make here) effectively restrict the generality of the treatment, therefore 
potentially introducing modeling inaccuracies. Both types of issues typically produce 
simpler expressions and algorithms than the ones discussed here, being probably this 
practical benefit the original reason behind them. The most radical example of this is 
the case of a freely jointed chain with constrained bond lengths and bond angles, and 
no other potential energy than the harmonic oscillators in charge of enforcing these 
constraints [SHf87|71l77ll33 134J . In such a situation, all terms in Vp(s) become inde- 
pendent from the conformation s (the torsional angles in this case), except for S* (s); 
and therefore all the difference between the rigid and stiff models can be attached to 
det g, which, in addition, can be very rapidly computed in terms of banded matrices 
[8W7] , 

Finally, it is worth remarking that, apart from the works that assume some of the 
approximations discussed above, the vast majority of constrained MD simulations in 
the literature do not consider any correcting term to Vs at all; sometimes giving 
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analytical or numerical reasons for it, sometimes simply not mentioning the issue 
[l23ll25l43l78l6lll26l62l58l4<)l59ll51ll3^l9ll5Dj . 



3 Numerical examples 

In this section, we provide a numerical example of the application of the general 
formalism introduced in this work to a real molecular system, using a classical force 
field. 

It is worth mentioning at this point that many numerical studies about the equi- 
librium of constrained systems exist in the literature, each one of them containing 
some of the approximations which we carefully enumerated in the previous section, 
and which are avoided in this work. 

Regarding the articles in which the stiff and rigid equilibrium distributions are 
compared, we can mention ref. |75j , where a simplified model of A^-butane is studied 
and the Fixman potential is found to be negligible; they use hard constraints, and 
they neglect both the Hessian and the s-dependence of detG, effectively using the 
(0,0)-stiff model introduced in the previous sections. In ref. [55], the same molecule is 
studied, under the same approximating assumptions; they find that the probability 
of the trans-gauche transition state is changed in a 20-30% by the Fixman correction. 
TV-butane is also analyzed in ref. |77j . where they do not see an appreciable effect in 
transition rates or relaxation times if the Fixman potential is included; they assume 
again that the Hessian and the s-dependence of det G are negligible, and they impose 
hard constraints. Under the same assumptions, in ref. |133j . the equilibrium of N- 
butane is shown to be altered by the inclusion of the Fixman correction if bond 
angles are constrained, but not if only bond lengths are (see below for more works 
discussing the bond angles issue). In ref. [73) . the equilibrium distribution of certain 
angles in simplified models of three and four beads is found to be different between 
the stiff and rigid models. The same four-beads model is simulated in ref. [7T], where 
the same comparison is performed and the rigid-plus-Fixman simulation is shown to 
be both equivalent to the unconstrained one, and different from the rigid simulation 
without correcting terms. In both ref. [73] and ref. [71], the simplicity of the model 
and the potential energy used to describe its behaviour make the hard and (0,0) 
assumptions exact. In ref. |86j . a linear freely jointed chain is considered, and the stiff 
torsional angles distribution is shown to be recovered from rigid simulations if the 
Fixman potential is included, which they find to have a non-negligible effect; they 
assume again the approximations that produce the (0,0)-stiff model, and they use 
hard constraints. The same system is studied in ref. [87], where Fixman's potential 
effects are measured as a function of the chain length, finding it non-negligible; again, 
the Hessian and the dependence of det G on s are not considered. 

A different but related family of works have numerically compared flexible and 
hard constraints. Among them, we can mention ref. [64] . were they measure appre- 
ciable discrepancies between the two types of constraints in the calculation of the free 
energy difference between liquid water and methanol. In ref. [49], the unconstrained 
dynamics is taken as the reference to which the flexible and hard rigid models are 
compared in the computation of the velocity autocorrelation function of liquid water; 
the authors find that the predictions of the unconstrained and flexible rigid models 
are very similar, and different from those of the hard rigid case. In ref. |50| the same 
conclusion is drawn. In ref. [41 , a simple four-beads system is studied and it is shown 
that the unconstrained average variation in bond angles is recovered if flexible con- 
straints are used. A similar study in A^-butane can be found in ref. |69j . where the 
unconstrained distribution of the central torsion angle is shown to be well reproduced 
by the flexible rigid (plus Fixman) simulation, but not by the hard one. Finally, 
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Fig. 4. The methanol molecule studied in this section. The only unconstrained coordinate, 
the dihedral angle (pe, is indicated with a blue arrow. 



in ref. |141) . flexible constraints are proved to better reproduce the low- frequency 
part of the vibrational spectrum than hard ones in a simple toy model. The reader 
should notice that only in ref. 69J the Fixman potential is used to correct the rigid 
dynamics, which makes the agreement between the flexible rigid simulations and the 
unconstrained ones in the rest of works somehow surprising, probably suggesting that 
the Fixman correction is not significantly important in the particular systems that 
have been explored. 

It is also worth to mention at this point the general agreement in the literature 
about bond angle constraints changing too much the equilibrium distribution with 
respect to the unconstrained case. Many works have studied this issue and have rec- 
ommended not to constrain bond angle coordinates on the basis of their numerical 
findings 76 58 59 90 132 91 . However, all the simulations in these works have been 
performed (i) using hard constraints, and (ii) without including the Fixman cor- 
rection. We have showed in this work that both these two assumptions should be 
considered approximations that might potentially compromise the accuracy of the 
final results. Therefore, we agree with the author of ref. [55] in suggesting the pos- 
sibility that, perhaps, constraining bond angles can be made accurately if flexible 
constraints are used and the most general form of the Fixman correction introduced 
in this work is used to produce stiff averages (the numerical cost of such an approach 
is very relevant too, but it is independent from the accuracy issue). The checking of 
this hypothesis in biologically relevant molecules is currently in progress and it shall 
be reported elsewhere. 

In the example calculation presented in this section, we have deliberately chosen 
a simple system in order to illustrate the theoretical concepts introduced in the pre- 
vious sections. We study here the methanol molecule schematically depicted in fig. 4j 
and we compute its potential energy using the AMBER 96 force field |152ll53j . The 
numeration of the atoms and the definition of the internal coordinates follow the 
SASMIC scheme, which is specially adapted to deal with constrained molecular sys- 
tems |45j . All the internal coordinates have been constrained, except for the principal 
dihedral angle, ipe, which then becomes the only variable parameterizing the inter- 
nal constrained subspace S. In order to produce the conformations in which all the 
quantities introduced in this work have been measured, we have scanned this dihedral 
angle from 0° to 180°, in steps of 10°. At each point, we have generated the hard con- 
formations simply by setting all the constrained coordinates to the constant values 
appearing in the force field, and the flexible conformations have been produced by 
minimizing the potential energy with respect to the constrained coordinates at fixed 
(pe- These minimizations in internal coordinates, as well as the computation of the 
potential energy Hessian, have been performed with Gaussian 03 |154j . All correcting 
terms to the different statistical mechanics models discussed in the previous section 
have been computed using home-made programs. 

First o f all, we measured the most accurate version of the flexible stiff free energy in 
sec. 2.4.3 the quantity F$ (s), as a function of (p^, as well as the Fixman potential, 
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Fig. 5. Free energy, _F S ' (s), of the flexible (2,2)-stiff model, and Fixman's potential, 
Vp(s), as a function of the dihedral angle </pg in methanol. Both quantities have been added 
an irrelevant energy reference in order for them to have zero average. 
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Std. a Max. 6 Std. a Max. 6 



Vs 
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0.37164847 
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-TS r k 


0.01116055 


0.03295965 


0.00000028 


0.00000098 


-TSl 


0.00008127 


0.00022596 


0.00000000 


0.00000000 


-TS C H 


0.01079596 


0.03198016 


0.01071757 


0.03179323 


-TS C A 


0.01078666 


0.03146866 


0.01070879 


0.03176306 


c/i 0) 


N/A 


N/A 


0.00358154 


0.01322298 


(7e (1) 


0.00000004 


0.00000013 


0.00358563 


0.01325454 


C4 (2) 


0.00000004 


0.00000013 


0.00357814 


0.01322534 



Table 2. "Standard deviation, and ^maximum variation in the conformational space of the 
methanol molecule of the different terms introduced in sec. |2.4| which affect the equilibrium 
free energy of constrained models. All values are in kcal/mol. 



Vp{s), in eq. (1261, which quantifies the difference between the stiff equilibrium and 
the rigid one. The results, depicted in fig. ml indicate that the variations in Vp(s) are 

(2 2) 

at least two orders of magnitude smaller than those in F s ' (s), and indeed much 
smaller than the thermal noise RT ~ 0.6 kcal/mol at room temperature, allowing us 
to conclude that, in this case, it is safe to neglect the Fixman correction. However, if 
we take a look at the variation of the different terms involved in the definition of Vf(s) 
in tab.[2j we see that —TS^(s) and —TS^(s) show a variation which is one order of 
magnitude larger than that of Vf(s). What is happening is that these two quantities 
are much correlated in methanol and have opposite signs, thus almost canceling out. 
Preliminary data in larger systems indicates that this feature is only present in small 
molecules such as the one studied here. Also, the significance of Vp(s) appears to 
grow with molecular size. Therefore, the particular results in methanol should not be 
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Fig. 6. Assessment of the hard rigid model for the methanol molecule, (a) Comparison 
between the flexible (2,2)-stiff free energy, F s ' 2 ' 2 ' (s), and the hard rigid one, F r (s) as a 
function of the dihedral angle ipe. (b) Difference between the two free energies in (a), as well 
as the difference between the potential energies, Vs(s) — Vs(s), and the correcting terms 
involved in the two models, (b) Difference between the flexible rigid free energy, F T (s), and 
the hard rigid one, F T (s). Again the difference between the potential energies, Vjj(s) — Vs(s), 
and the difference between the respective correcting terms are also depicted. 



extrapolated in face value, but rather seen as an illustration of the concepts involved 
(probably the same can be said about the works discussed at the beginning of this 
section which analyze small toy systems and TV-butane). 

The difference between flexible and hard constraints in this molecule is larger than 
the difference between the stiff and rigid models, as we can appreciate in fig. [6j where 
the hard rigid model is compared both to the flexible (2,2)-stiff one (fig. [6]d) and to 
the flexible rigid one (fig. |6j:). In both cases, the maximum variation of the difference 
in the conformational space of methanol is approximately 0.04 kcal/mol, an order of 
magnitude larger than the variation of Vp(s) in fig. [5] We can also see that both the 
difference between the flexible and hard potential energies, Vs(s) — Vs(s), and the 
difference between the correcting terms in the respective models contribute to the 
total discrepancy between the flexible and hard cases, being the second contribution 
somewhat larger than the first one. The fact that the graphs in figs. [6]d and [6]: seem 
identical to the eye (they are not exactly the same numerically) is again caused 
by the already mentioned correlation between the terms —TS^(s) and —TS^(s) in 
methanol. Preliminary data in larger molecules suggest that these two comparisons are 
different in general, and also that the difference Vs(s) — Vs(s) becomes increasingly 
more important with system size if the bond angles remain among the constrained 
coordinates. This is probably due to the possibility of more steric clashes in longer, 
more flexible chains, and it is in agreement with the works discussing the bond angles 
issue that we have summarized at the beginning of this section. 
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Fig. 7. Variation of the constrained coordinates in the flexible case in methanol, as a 
function of the unconstrained dihedral angle ip&. (a) Variation of the bond lengths, (b) 
bond angles, and (c) phase dihedral angles with respect to the constant values appearing in 
the force field and defining the hard models. 



Finally, we can see in fig.[7]that the constrained coordinates depend on the confor- 
mation in the more general flexible expected from the discussion in sec. |2.3.2| 

4 Conclusions and future lines of research 

In this review, we have attempted to provide a unifying view of the statistical equi- 
librium of constrained classical mechanics models of molecular systems. To this end, 
we have introduced the most general formalism (i.e., including the fewest approxima- 
tions) compatible with the possibility of having analytical results. From this advan- 
tageous standpoint, we have been able to rationalize most of the previous works in 
the literature, clearly identifying the underlying (often implicit) assumptions behind 
the different analyses. Also, we have tried to provide the reader with a coherent vo- 
cabulary which we hope facilitates future comparisons and reviews. Finally, we have 
shown a practical example of the different theoretical concepts. 

Along the text, we have suggested a number of possible future lines of research that 
we believe could be interesting for distinct reasons, and some of which we have already 
started to pursue. We conclude this account by briefly collecting and commenting 
them: 

— Define the constrained coordinates in an adaptive way, instead of simply picking 
whole sets of internal molecular coordinates of the same type (e.g., all bond angles 
involving hydrogens). If properly done, such an scheme could not only help to 
make calculations more accurate, but also to assess how good are the simpler 
approaches followed at the moment in the literature. 
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Develop a rigorous and general quantum mechanical derivation of constrained 
molecular models that could help decide which one of the classical models dis- 
cussed in this work is more appropriately justified on physical grounds (if any of 
them) . 

Compute the equilibrium PDF of the non-Hamiltonian, flexibly constrained dy- 
namics described in ref. [13], maybe using the technique in ref. (115] . 
Compare the option chosen in this work for defining the flexible constraints based 
on t he minimization of the potential energy V(w) to the alternative discussed in 
sec. 2.3.2 which uses the free energy F(q) that appears in the exponent of the 
marginal PDF of the unconstrained dynamics in the space of the positions. 
Use techniques for the computation of free energy differences in order to numeri- 
cally calculate the marginal PDF of the unconstrained dynamics in the space of the 
unconstrained coordinates u [see eq. (85)], and compare it to the approximation 
to it provided by the stiff model. 

Solve the discrepancy between the correction to the rigid dynamics found from 
statistical analyses like the one in this work (i.e., the Fixman potential), and those 
obtained in some studies based on dynamical considerations, such as ref. [30] • See 
also sec. 12.51 for further details. 

Extend the analysis in sec. [3] to larger, more relevant biological molecules. 
Explore the hypothesis, mentioned in sec.[3j that flexibly constraining bond angles 
and including the most general form of the Fixman potential to correct the rigid 
equilibrium may result in more accurate constrained simulations than the ones 
reported so far in the literature, where constraining bond angles is typically not 
recommended. 
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